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ON DIFFUSION BY DISCONTINUOUS MOVEMENTS, 
AND ON THE TELEGRAPH EQUATION 
By 8. GOLDSTEIN 
(Institute of Technology, Haifa, Israel) 
[Received 8 August 1950] 


SUMMARY 

At time ¢ 0 a large number of non-interacting particles start from an origin and 
move with a uniform velocity v along a straight line for an interval of time r. To begin 
with, half move in each direction. Thereafter, and at the end of each successive 
interval of time 7, each particle starts a new partial path; it still moves with speed v, 
and there is a probability p that it will continue to move in the same direction as in 
its previous path, and a probability q ( 1 —p) that the direction of its velocity will 
be reversed, so the directions in any two consecutive intervals are correlated with a 
correlation coefficient ec = p—g. The partial correlations for non-consecutive inter- 
vals are zero. The difference equation is found for the fraction y(n, v) of the number 
of particles at a distance y vor from the origin after a time ¢t nt; it is shown how 
y(n,v) may be computed; asymptotic formulae for large n are found, both for a fixed 
value of p/q and for a fixed value of ng/p. The limiting density distribution (and the 
limiting characteristic function) are found when n — 0, tT —> 0, with nr t, vuT Y; 
and in the limiting operation c |—7/A, with A constant, so that c — 1, and the 
speed v is kept constant; the limiting form of the difference equation is the telegraph 
equation, with » = (LC)-!, A L/R, where L, C, and R are the self-inductance, 
capacitance, and resistance per unit length; and the limiting density distribution is 
the solution of this equation for an instantaneous source. If p+q + 1, and there is, 
at the end of each interval of time 7, a non-zero probability, |1—p—q, that a particle 
will escape from the system, then the same limiting operation, with 7/(1—p—q) = G, 
G constant, leads to the telegraph equation with leakage, the leakage resistance 
being G/C. The solution of the telegraph equation is further considered (and in 
particular is given in terms of Lommel’s function of two variables) when one end of 
along cable is held at a constant potential for t > 0. 


1, Introduction 

BEFORE considering diffusion by continuous movements, G. I. Taylor (1) 
ina paper published in 1922, considered briefly an extension of the usual 
treatment of discontinuous random migration in one dimension. The 
problem is phrased as a physical problem in diffusion by considering that 
at time ¢ = 0 a large number of non-interacting particles start from an 
origin and move with a uniform velocity v along a straight line for an 
interval of time +, during which each moves a distance d = vr; to begin 
with, half move to the right (in the positive direction along the line) and 
half to the left (in the negative direction); after time 7, and after each 


interval of time 7 thereafter, each particle starts a new partial path, and 


again moves with velocity v through a distance d in time 7, but may either 
continue in the same direction as in the previous path or may move in 
[Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 2 (1951)] 
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the opposite direction. The process continues for a large number of partial 
paths and intervals 7. In the usual theory of random flights there is no 
correlation between the directions of motion of a particle in any two time- 
intervals. If, however, the direction in any one interval [mr, (m-+-1)r] is 
correlated with that in the next interval {(m-+1)r, (m+ 2)r] with a correla- 
tion coefficient c, and the partial correlations for non-consecutive intervals 
are all zero (so that the correlation coefficient between the directions in the 
intervals [mz,(m+1)r] and [(m-+s)r,(m+s-+1)r] is c’), then Taylort 
showed that the mean square distance from the origin after time nz = T is 

yaftte pp 2e(1—e") 9) (1) 

\1—c (l1—c)® J 

As 7 > 0 and n> ©, with v fixed, the diffusion approximates to a con- 
tinuous diffusion of some kind; the mean square distance stays finite if 
7/(1—c) > a non-zero limit. If +/(1—c) = A, then in the limit the correla- 
tion coefficient between the directions of the velocity of a particle at two 
times separated by an interval t is e~’4, and the mean square distance from 
the origin after a time 7’ is 

v{2A T—2A2(1—e-74)}, (2) 
which becomes v?7? for small 7’, as it should, the correlation not having 
fallen from unity. 

These are Taylor’s results. Taylor went on to consider diffusion by 
continuous movements; the process considered above cannot lead to 
turbulent diffusion. Some years ago (1938/9) I considered more fully the 
density distribution of the particles after n steps, and the approach to the 
limit. This work is reported here. (It appears from my notes that at that 
time Mr. J. Wishart had obtained independently the formula (16) below 
for the moment-generating function, and had drawn my attention to a 
paper (2) he had published on a similar problem.) 

After the difference equation for the fraction, y(n,v), of the number of 
particles at a distance vd from the origin after a time nz, and the formula 
for the moment-generating function, have been found in section 3, the 
computation of y(n,v) is considered in sections 4, 5, 6, where asymptotic 
expressions are found for large n, both for a fixed value of ¢ (which may 
be positive or negative) and for a fixed positive value of n/r, where 
r = (1+c)/(1—c). The passage to a continuous distribution, when n > ©, 
+>0, d+0, c>1, with v constant, nr =t, vd = y, 7/(1—c) =A, 
n/r > B = t/2A, A constant, is considered in sections 7 and 8. The limiting 
density distribution, and the limiting characteristic function, and their 


+ c®? occurs in place of c in the numerator of the second term in ref. 1, It appears that 
on the last line of p. 199 of ref. 1, n should be replaced by n—1. 
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connexion, are considered in section 7; in section 8 it is shown that the 
limiting density distribution satisfies the telegraph equation (without 
leakage), and with v? (LC)-1, A L/R, where L, C, and R are the 
self-inductance, capacitance, and resistance per unit length; this equation 
is the limiting form of the difference equation for y(n, v); and the limiting 
density distribution is the solution of this equation for an instantaneous 
plane source at the origin. The conditions under which the limiting opera- 
tion is carried out should be stressed: the partial correlations for non- 
consecutive intervals are all zero, c > 1 in a fixed manner (c = 1—7/A), 
und the particle speed is constant. (With three possible velocities, for 
example, +v and 0, a third-order equation is obtained (section 10).) In 
section 9 it is shown how a solution of the telegraph equation when one 
end of a long cable is suddenly raised to, and thereafter held at, a constant 
potential is related to the solution for an instantaneous source, and how 
various forms of the solution of that problem are related among them- 
selves, with certain mathematical relations involving Bessel functions and 
Lommel functions of two variables. In section 10 some further develop- 
ments are mentioned; in particular, it is shown that if the sum, p+q, of 
the probability p that a particle will continue to move in the same direction 
and the probability ¢ that it will reverse its direction of motion, is not equal 
to unity, but there is a non-zero probability, 1—p—gq, that it will escape 
from the system altogether, then the limiting form of the difference equa- 
tion when n > 00, tr > 0, with nt = t, 7/(1—c) = A, 7/(l—p—q) = G, is 
the telegraph equation with leakage, the leakage resistance being G/C; the 
density distribution function is easily found for this case. The solution for 
onstant potential at one end again involves Lommel’s function of two 
variables. 


2. Notation 


Ay, A, for definition, see equation (18). 
b,, by for definition, see equation (48). 
be for definition, see equation (52). 
¢(= p—gq) the correlation coefficient between the directions of the 
velocities in two consecutive partial paths. 
d (= vr) the length of each partial path. 
f f(x) { eiuv dF = M(n,ud). 
y=—a 
ins the limit off for the ‘continuous’ distribution (equation (94)). 
ky, ke for definition, see equation (19). 
LL, 1, for definition, see equation (47). 


for definition, see equation (52). 
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m an integer. 

n the number of partial paths. 

p the probability that a particle continues to move in the 
same direction at the end of a partial path: the Heaviside _ | 
operator in section 9. 

q the probability that the direction of motion of a particle is | 
reversed at the end of a partial path. | 

r= ply 

8 an integer. 

t time; used in place of 7’ for nz from section 8 onwards. 

u a real variable. 

Uy sin-!(1/r). 

v the uniform speed in each particle path. 

Ww (nr)tu. 

y (= vd) distance from the origin. 

A 7/(1—c) (when 7 > 0,c > 1). 

A, defined in equation (66). 

B=n/1 

CU capacitance per unit length. 

D nd VT’. 

F F(y) is the probability that the coordinate of a particle 
relative to the origin is < y. 

f* the limit of F in the continuous distribution (see equa- 
tion (95)). 

G = 7/(1—p—q) when there is leakage (section 10). 

L,, i, Bessel functions of imaginary argument (equation (81)). 

K a constant (= v/n in (63) and in section 6). 

kK, a constant. 

L for the meaning of L in equation (66), see equation (65). 
In sections 8, 9, 10, L is self-inductance per unit length. 

n 
M M(n,u) is the moment-generating function }  y(n.v)e"". 
yn 

N half the total number of particles (section 8). 

R resistance per unit length. 

S 1/S is the leakage resistance per unit length. 

S,, Sy, Sg, Sy, S; integrals defined in equations (42), (43), (44), (54) 
and (60). 

7, = oF (from section 8 onwards, ¢ is used in place of 7’.) 

U Lommel’s function of two variables (ref. 5, sections 16.5 


16.59, pp. 537-50). 
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W (vt—y)/(2Av). 
9 r | | | - Uj 
th Ms Me ay “G (! } 
side y B(1—v?/n?) B(1— K?)! Bi —y?/D?)! (= (vt?—y?)!/(2Av) in 


the limit). In section 10, Y is the positive square root of 


. (1 (e2) 
4\A G v- 

(n,v). B(n.v) the fraction of the number of particles at a distance vd 
from the origin after a time nr, and arriving from the 
left (a(n,v)) and from the right (8(n,v)), respectively. 

y(n, v) (= a(n,v)+A(n,v)): the total fraction of the number of 
particles at a distance vd from the origin after a time nr. 


4 for definition, see equation (55). 
v temperature. 

K thermometric conductivity. 

v an integer. 


p l/r q/p- 
line density. 


T time interval to describe a partial path. 


db for definition, see equation (56). 
rucle x y Dsin x, Y B cos x (equation (91)). 
us sin-l(r sin w). 

qua 
3. The difference equation for the density of the distribution. 

" The solution for the moment-generating function 

Let p be the probability that a particle will continue to move in the same 

direction at the end of a partial path, and q = 1—p the probability that 

Rx the direction of its velocity will be reversed. The correlation coefficient is 

neti Cc Pp q- (3) 
' . . . . . 

)e Let y(n,v) be the fraction of the number of particles at a distance vd from 
the origin after a time nr, a(n,v) being the fraction that, moving towards 
the right (in the positive direction), arrives from the left, and B(n,v) the 
fraction that arrives from the right, moving towards the left. (n is positive; 

n4 vmay be positive, negative, or zero.) Then, by elementary reasoning, 
' 
y(n, v) x(n, v)+B(n, v), (4) 

16.5 x(n, —v) 3(n.v), y(n, —v) y(n, v), (5) 

ro 
x(7, v) 8(n,v) = y(n,1 0 when n—v is odd, when vp - n, 


and when vp 
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x(1,1) = B(1, —1) = 3, x(1, —1) = B(1,1) = 0 
y(1,1) = y(1,—1) = 3 ° 
x(2, 2) = B(2, —2) = dp, x(2,0) = B(2,0) = 4 
x(2, —2) = B(2,2) = 0, : (8) 


y(2,2) = y(2, —2) 
and the following recurrence relations hold: 
x(n+1,v) = pa(n, v—1)+qB(n, v—1) py(n,v—1)—cB(n,v—1), (9) 
B(n+-1,v) = pB(n,v+1)+qa(n,v+1) = py(n,v+1)—ea(n,v+1). (10) 
From (9) and (10), 
x(n, v-+1)+B(n, v—1) pa(n—1,v)+qB(n—1, v)+pB(n—1, v)+qa(n—1, v) 
y(n—1,v). (11) 
Hence, by addition of (9) and (10), 
y(n+1,v) = ply(n,v—1)+y(n,v+1)]—cy(n—1, v). (12) 


This is the difference equation for y(n,v). Now define 
n . n 
M(n, u) > y(n ve" = y(n, 0)+2 > y(n, v)cos vu, (13) 
v n v=1 
where wu is a real variable. From (7) and (8), 
M(1,u) = cosu, M(2,u) = q+pcos2u 1—2p sin?w. (14) 


Multiply (12) by e’’” and sum for all values of v from —(n+1) to n+. 
Since y(n,v) = 0 for v - n and for v > n, the result is 


Vv n 


1 
M(n -_ B wu) pe" > y(n, v- L)et Du_ 
1 


Vv n 
_ v+tl=n : poRn—s : 
+ pe iu > y(n,v | 1 )e i(v+lu_e > y(n 1, v)ei”™ 
v+1 n v (n—1) 
= p(e'"+e-) M(n, u)—cM(n—1, u) 
2p cosuM(n,u)—cM(n—1, u). (15) 
This is the difference equation for M(n,u). Its solution is 
M(n,u) = k,a®+k, a, (16) 
where a, and a, are the roots of the quadratic equation 
x?— 2px cosu+e = 0, 
: ‘ ) r ] r—l S 
i.e. with r= I ; p= ‘ q - c= , (17) 
q r+] r+] r+] 
a, = peosu+q{1—r? sin2u}}, ad, = peosu—q{1—r*sin?u}*. (18) 


(This result may also be obtained by manipulation of the transition matrix 
of the Markoff chain: ef. ref. 3.) 
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The values of k, and k, are found from those of M(1,u) and M(2,u) in 
(14), and are 


l -O8 ] :08 
ky 14 one ; els are wi | (19) 


2} © (1—?r?sin2w)! a r? sin?u)! 


- 


The moments of the distribution about the origin are found from the 
expansion of M(n,u) in ascending powers of u; the odd moments are zero, 


n 


and > v?"y(n,v) is (—1)™ times the coefficient of u?"/(2m)! in the ex- 
yon 
pansion of M(n,u). For example, 
> v*y(n, v) nr—4(r?—1)(1—c") (20) 
(in agreement with (1)) and 
> v4y(n, v) 3n?r?— 6nr3 + 4nr+ 3(r?— 1)(9r?—1)(1—c") —3nr(r?—1)e”. 


(21) 
If we carry out the limiting operation +t > 0 with nt = T, 7/(l—c) = A, 
as before, then rr > 2A, and we have, as before, for the mean square 
distance from the origin, 
d? > v*y(n, v)—> v>| 2A ‘i -2A( l—e it ) | (22) 
also 
d* ¥ v4y(n,v) > 12v4{ 42722437 (2+e-7/4) 4 6A4(1—e-7/4)], (23) 
and we verify that this becomes v‘7" for small 7’, as it should. 
It follows from (13) that, in terms of M(n), 


I; 
y(n, v) M(n, u)cos vu du. (24) 
. 0 
since 
a,(7—u) a,(u), k,(7—u) k,(u), ! (25) 
20 
’ 
M(n.7—u) (—1)"M(n, u) 
y(n.v) = 0 if n—v is odd, as it should, and if n—v is even 
T 7 
9 oe FF 
y(n, v) M(n,u)cosvu du = | k, af cos vu du. (26) 
7 * 7 * 
0 0 


4. The computation of y»(n,v). Asymptotic expansions for n large, 
r fixed 
The values of y for n | and 2 are given in (7) and (8), and the recur- 
rence relation in (12); for small values of n it is easy to write down the 
values of y explicitly. For any value of n, 


y(n, n) ig*, (2 


) 
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and for an integral s > 0 and < 3n 
n—3(> l ] x a re 
atthe — 1) i. holt 
’ (r+1)"-! ane 2! 2! r4 
(s—1), (n—s—1), 1 
3! 3! 76 
n—-3 ] ] re 4 ra a 
nt it~ tia ~o~¥h— ( 1) (n 1), | 
Ar+1yap 7? Pe'3 2 21 orl 
(s— Ta ( a Hs | oofy (20) 
i 8! 5 ‘ia 
where (si. x(x—1)...(2—m-+1). (29) 


If s is not large, the values of y are easily found from this formula for any n, 
When r = 1 the problem is the classical one of random flights in one 
dimension; the solution is 
l n! 


on Slim —v\t! 51 ot? (90) 
2” 13(n—v) shy 9(m+v)3! 


y(n, v) 


with the asymptotic expansion (4) for large n and v? = O(n) or o(n) 


9\1 

= \5 ,2/6 
y(n, v) pe =e—v*/2n vy 

7H 


1 ; “( es ] ") ] (1 _44 y2 ‘ 38 vw 12 : 1 v a 
4n wry 


n 3n?] ° 32n? 3 n 3 nn? 5 n?' 9n! 


When r + | an asymptotic expansion similar to (31) may be found from 
the integral (26). ((28) may be expressed in terms of hypergeometric func- 


tions, but it seems simpler to use (26).) The result is 


f 2 \3 , 1 /r?+1 v= 3r?—1 v4 
y(n,v) ~ * e-vilaen| y Or 4 
7rn 4rn\ 2 n 6r? on? 


l 2 10727+-3 1874+ 30r?—4 vp? 8174— 107?-+5 v4 


327727 \ 4 or n 6r? n* 
« 
Mr S41 @ sf 2.2 2 
« 207r4—10r?-+-2 v8 (372?—1)? | (32) 
5r3 n3 3674 on! 


A comparison of numerical values will be set out before a brief descrip- 
tion of the derivation of the formula. 

For numerical comparison, » is taken as 15, and r as 4, 1, 2, and 5. In 
Table 1, each column headed I contains the correct values of y(15,v); the 
column headed II contains 107 times the error in the asymptotic expansion 
(32) taken as far as the terms shown there explicitly—i.e. 10° times 
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column | minus values calculated from the terms shown explicitly on 


») 


the right-hand side of (32)). 


TABLE | 





u <,¢ i ? 5, 3 

I] | I] I I] | II 
25 0°1405597 2¢ 0°05744006; 625 
| 77 527405 17 O'1 249551 d5O0 0°0542032 1161 
} 14 72 0°098 3322 170 | 0°0779941 2026 
” 52 0°00760935 ISs9 0°069 3103 2974 
69 "O2QQS551 788 0°0588573 2482 
47 070194525 +613 0°0474173 10g! 
j 15 35 0°007 2792 554 0°0355275 + 7095 
14 0°0017127 1311 0°0 359433 + 155391 


For the values shown in the table the accuracy of (32) decreases, for a 


fixed n, as 7 increases from The formula is not satisfactory, for any n, 

r is too small or too large. The case of small r is not further considered 
here. For large r and large n the present method would not normally be 
ised (see section 6 and Table 2). 

The derivation of (32) from (26) is fairly straightforward, but a satis- 
factory estimation of the error terms is rather tedious, since the methods 
to be used and, to some extent, the results, depend on the value of r. Care- 
ful estimations of the error terms appear to be of no great interest, and are, 


for the most part, omitted; only a brief sketch of the proof of (32) will be 


given. 

Let us begin with the case 7 1, and use the second expression in (26). 
For r |, a, is real and positive, and decreases from | to —c as u increases 
from 0 to 7, where —c (l1—r) (1-+-r); k, is positive and decreases from 1 


0, 80 hk, a} also is positive and decreases from | to 0. 
We may, for a certain range of values of u, expand loga,+4ru? in 


powers of u 
loga, 5rue r(3r2— 1)u4 aor (45r4 30r?+ 1)u®+ O(u8). 


We may find an expression for a e!"""’ by expanding exp{n(log a,+ 4ru?)} 
n powers of n: if error terms are included, we may then substitute for 
ga,+ ru*. To include satisfactory error terms we require that 


exp{n(loga, 4 tru*)} 


should be bounded, so we restrict uw to the range 0 < u < K,n~', where 


K, is a constant; that exp{n(loga,+4ru?)} is then bounded follows from 


5 


d rsin u ; rue ' 
log a,) a rsin u - ra +-—_, (33) 
du (l1—r? sin?u)! 6 
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"4 
o rut _ rKj 


so loga,+4ru< e aa, 34 
—— 24 ~ 24n -_ 
We may now use finite expansions with error terms to show that, if 
w = (nr)tu, (35) 
then 
2 A A 
i we. we u 
k, ate L - (r2— 1) — - (3r2— 1) + —_. (r?—1)(9r?—1)— 
4nr 24nr 48n*r* 
ws or 9 ws 9.9 9 
— —. (135r4— 120r?+ 1) + ————. (3r?—1)?+..., (36) 
1440n*r* 1152n*r 


with errors of orders w®/n?, w’/n3, w!9/n3, w!?/n3, and with w restricted to 


the range 0 < w < K,r'n'. 
The whole of the asymptotic expansion comes from 
K,in-t 
a, 
k, a} cos vu du; (37) 
7 


a 


0 


k, < 1 and a, is a decreasing function of uw, so in the rest of the range of 


integration the integrand has a modulus less than the value of a? at 
u = K,n-, which in turn is less than exp(rK}/24)exp(— 4rK?n!); and the 
error due to neglecting the rest of the integral is less than twice this. 

In (37) we now substitute from (36), and change the variable of integra- 
tion to w. The range of integration is then from 0 to K,r'n*. The range 
is extended to be from 0 to 00, the error so introduced being of order 
n-texp(—4rKjn') at most. The value of k,a? from (36), multiplied by 
cos|vw/(nr)*], is integrated term by term, the equality 


Oo o 
‘ad e ° 9 
se vw . * yw T ae 
e hu “yp2m cos . dw e—}v?/nr re | e—}2 x + ; dx (38) 
(nr)? 7 (nr)? 





ry 


0 0 

being used. The error in stopping after a finite number of terms is of the 
order of the first neglected term. Thus (32) is obtained for r < 1. 

For r = 1, k, is 1 in (0,47) and 0 in (47,7); a, is cosu in (0,47) and 
0 in (47,7). Equation (30) and Stirling’s formula lead at once to the 
result (31), but the calculations from the integral could proceed as before: 
however, in the expansion of exp{n(log a,+ }u?)!, wu need not be restricted 
to an interval (0, K,n-'), but may be allowed to have any value in an 
interval (0, 47—8) for any positive 5, since loga,+ 4u? is negative, and the 
exponential therefore bounded, in the whole of that range. The error 
terms are correspondingly altered. [A somewhat similar remark applies to 
a certain range of values of r < 1; if r is not too small, there is an interval 
of u, independent of n, in which log a, + }ru? is negative, which may replace 
the interval (0, K,n-*).| 
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, Consider now the case r 1, which is the case of greatest interest. Let 
(34) 
cam 
U, sin , (39) 
> 
(35 Then k, and k, are infinite at vw = u,; k, increases from 1, and k, decreases 
from 0, at wu = 0; a, and a, are finite, and each is equal to c!, where 
Cc (r—1)/(r+1), 
at u u,; @, decreases from 1, and a, increases from c, at u 0, to this 
(36 ee ji ? 
value. M(n, uw) is finite and equal to c!"(1+n/r) at u u,; it decreases to 
d to this value from 1 at wu = 0. 
Between wu, and 37, a, and a, are conjugate complex, and so are k, and 
k,, with 
l ee ] 2 cos u 
a, [7 cos u+i(r? sin?u— 1)! |, k, 1—i —~ . 
(37) rt] 2 (r* sin?’u— 1)? 

(40) 
ye of M(n,u) is oscillatory, and it may be proved by induction that M(2n,u), 
* at M(2n+1,u)/(cosu), are both polynomials of degree in cos*u, each with 
1 the n zeros between u u, and u Ar. 

From (26) dary(n, v) S,4 Sy +S, (41) 
ora 
where " 
ange 
y ° Na = ») 
vailien S | k, a} cosvu du, (42) 
0 
d by 
; u 
S, | kya} cos vu du, (43) 
0 
(38 
und Sy | (k,af+k,az)cos vu du. (44) 
f the Write 
p : (45) 
) and r 
» the ub = sin-(rsinu), (46) 
fore 
| l COS us l cos os m 
icted l, : 1+ ae ‘ l, l Tr Foe |’ (47) 
2 (1— p* sins)! 2 (1—p? sin?ys)? | 
in an 
d the } | . ] 2 cjn2,/,)4 ] } I . 2 cin2,/,\4 
} 0, | p cos #+-(1— p* sinh)? ], Ne [ p cos s—(1— p* sin*y)! | 
erro! 1+ P Fi ] T p 
ies ti (48) 
; | > , . 
erval | (80 that /,, /,, b,, b,, are the same functions of p and % as k,, ky, a,, ag are 
place of r and uw); then 


k, du pl, dip, k, du pl, dip, a, b,. as -by, (49) 
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and - 
S; p 1, bt cos{v sin-!(p sin x)} dy, (50) 
0 
in 
Sy (—1)"+1p 1, bf cos{v sin~1(p sin x)} des 
0 
p | 1, b{ cos{vsin-(psin w)} dip. (51) 


Le 
Since as y& increases from $7 to 7, 1, b} decreases from its value 
1 
[(1—p)/(1+p)}!" 
at 4a to zero at 7, \S,| < $ap[(1—p)/(1+p)|*". We may now write a 
transformation of the integral, which we shall use later. If 


] sin ws } (1 — p? cos*ys)! — p sin ob (52 
. - ‘ de ; 52 
‘ (1— p* cosy)? ' (1—p?)! 
‘ l ae 
then Sy bo 4) Sy. (53) 
l+p 
iz 
where S, | 1, b& cos{v sin-!(p cos w)} dab. (54) 
0 


S, may be omitted in establishing the asymptotic expansion (32). To con- 
sider S;, commence by putting 


, (7? sin?u— 1)! 


@ = tan (= arga,), (55) 
reosu 
, COS Uu , ] ar 
d@ = tan ——- . = tan cot é} | arg k,). (56) 
(7* sin*u— 1)? r 
Then, from (40), 
sar 
; y—i1\i" i sin u = 
S3 : (r2— 1)! ——~  cos(n6—¢)cos vu du. (57) 
r+1 (7* sin? — 1)? 
Uy 
Now (r2— 1)cos?0 r2 cos?u, (r2— 1)sin26 r? sin?’u— 1, (58) 


so if the variable of integration is changed to 8, 


. 1—p\'" ~~ r 
™ | (1— p*)'S5, (99) 
l+p 
where se 
S; | cos{né—tan~1(p cot @)}cos{v cos (1 — p?)! cos 6} dé. (60) 
0 
S;| < 4, soS; may be omitted. The asymptotic expansion (32) therefore 


comes entirely from S,. More exactly, we shall prove later that 


28,+S8,=0 forv<n | 


' (61) 
in(l1—p)"*! forv=n>0] 
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SO _ 
lery(n,v) = S,—S, p | 1, b} cosfyvsin-“(psinys)} db forv <n : 
, ; (62) 
0 
S, Sy pa(1—p)” 1 for v n > VU. 
The asymptotic expansion (32) may be established fro 52) for r > 
ptot | ») may | tablished from (62) f l 
in the same way as it was established from (26) for r < 1, the only new 


calculation required being that cos~'{vsin-1(psin #)}, after the substitution 
» — (nr) for the variable of integration (corresponding to the substitu- 
tion (35)), must be expanded in powers of n-!. Alternatively, and rather 
more easily, the argument may start from the original integral (24), and 
proceed exactly as before; uw may be taken over the interval (0, u,—68), for 


Lru? is negative; and the transformations 


5 


any positive 6, since loga, 
shown above may be used to find an upper bound to the error in neglecting 
all but the integral, from 0 to u,—8, of the integrand in (42). 


The very large jump in the error in Table 1 for r = 5, between v = 13 
and v 15, is largely accounted .for by the term har(1 p)"+t in (62). 
The asymptotic expansion (32) is valid for v* O(n) or o(n); when v is 


nearly equal to n, direct computation from (28) is not difficult; when v/n 
is equal to a constant K, not nearly equal to 1, then from (30), for r = 1, 
2\} 1— K\'Kn 3+K2 1 
y(n.v) ~ | ime K*)- »| ; I == rae (63) 
nt 1+ kA | 12(1— K?) n | 
and for ? | this formula may be generalized, by use of the method of 


steepest descents on the integrals in (26) and (62), to 





. (1— Kk?) ]A"Y A, 
(l—A%) “118 ; oe (64 
L+k | n- : ) 
where L (K24 72(1— K2)i! (65) 
and 
, _ 37°(3r?—1)— K*(18r4— 18r2— 2) + K4(9r4— 15? + 6) 
i 24(1—K2) 13 
de il = (66) 
2 L*( L T r) 
5. Proof of the relation between the integrals S, and S, 
To prove (61), with p 1)? |, it follows from (53) and (59) that it is 
necessary and sufficient to prove that 
NS P N () tor pv 
*~ (1— pyr sialic . (67) 


l 2 


re | p)(1 p?)iu ) for v n mb) 
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S; is given by (60), and it is easily seen that 
47 


S. 4 re e in® dane 


o 


sin 6—ip cos @ " — ‘ : 
Ge a cae (| (1 —p?)? cos 6+ t(p? cos*6+ sin?) + 


+[|(1—p?)! cos 6—i(p? cos?6+-sin?6)*}"} d0. (68) 


S, is given by (54). In the integral, substitute 


(1—p?)! sinha, = psiny, (1—p2)!cosha, = (1—p2cos%p). (69) 
Then 
b, = e-%, a (pes sinh = 
(—9F pcosh x,+sinh x, , (70) 


exp|isin-'(p cos ys)] = (1—p?) cosh x, +-7i(p? cosh?x, —sinh2z,)+ 
and as % increases from 0 to $7, x, increases from 0 to tanh~'p, so 


tanh 1p 


vats) -_ a , 1 
- Sy = | é ns,(200e 7 —sinh “1 *{{(1—p?)! cosh 2,4 
(I1—p*) 2 F pcosh x,+sinh 2, 

) 


++ i(p? cosh*x, —sinh*x)* |»+-[ (1 — p?)! cosh x, —i(p? cosh*x, —sinh?x,)#]"} da. 





(71) 
To prove (67) consider 
1 f .. peoshz+sinh z : . 
: enz_P ——-~ {[ (1 —p?)! cosh z+-i(p? cosh?z— sinh2z)*]"+ 
2. (p? cosh?z— sinh?z)! 


+-[(1—p?)! cosh z—i(p? cosh2z—sinh?z)?]"} dz, (72) 


where z = x+iy, round the contour shown in Fig. 1. AB is x=0, 


0<y<}a; BC is y=}n, 0>2x > —X, where X >~; CD is at 
«x = —X; AED is along the negative real axis, and is indented at E, by 
a small semicircle whose radius > 0, where EF is at x = —tanh-'p. Then 

_; omemens —<_______————_B 

Y A 
cee : . { 
E 
Fic. | 


the real part of the integral along AB is S;; when, as here, $(n—yv) is an 
integer, the integral along BC is purely imaginary; the integral along 
CD +0 for v <n and > —}n(1—p)(1—p?)"-» for v = n; the integral 
along DE is purely imaginary; the integral round the small semicircle at 
E +0; and the integral along ZA is —pS,/(1—p?)#, according to (71). 
Hence (67), and (61), are proved. 
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6. The asymptotic expression for large n, n/r fixed 
Consider next the values of y(n, v) for large values of n and r, when r/n 
is constant, and equal to 1/B, say. Thus as noo, r>o,c—>1. Also 
let vn be constant, and equal to K; then K is the ratio of the distance 
18) y — vd from the origin to the greatest distance D = v7’, where v is the 
particle velocity and 7' the total time of motion, equal to nz or nd/v. 


Note also that to accord with Taylor’s notation, in which 7/(l—c) = A, 
)9) 


we must have T 
B , (73) 
A(1-+c) 

“()) As n> 00, we may suppose that 7 0 and d-> 0 in such a way that 
T = nr and D = nd remain finite. Then in the limit the mean square 
distance of the particles from the origin is given by 

y* | ] oR = 
o ,(1—e-™). (74) 
py? JT B 2 B? 
Cf. (2) and (22).) y(m,v) is given by (62), p l/r Bn, and as n > «x 

l, > 4(1+ cosy), b} > exp(— B+ Beosy), | (75) 

Uy. : Re ‘ . io 

a cosfv sin-}(psin¢s)! > ecos(K Bsin ws J 

1) t f f 
The integrand tends uniformly to its limit and 
nS, > 4Be-® | (1+ cosa&)exp( B cos ys)cos(KB sin p) db. (76) 
0 
(72) Similarly, _ 
0. nS, > 4Be-® | (1—cosis)exp(— B cos )cos(KB sin p) dif. (77) 
at 0 
; Hence, for v n 
by = 
hen brny(n,v) > Be-® | {cosh( B cos ys)+ cos % sinh( B cos %)}cos( AB sin) db. 
i (78) 
Sincet si 
») fe 
- cos( KB sin ys)cosh(B cos ps) db = 1,(Y), (79) 
0 
2 . 
cos(AB sin y)sinh( B cos )cos yy db = BL(Y)/Y, (80) 

is an . 

Lone where J), J, are the ‘Bessel functions of imaginary argument’, 

ooral | —. (12)2” = 4z)2m+ 

I,(z) S iam [,(2) Zz _\#) “Ft? (81) 
le at LL, (m!)? m! (m-+1) 

(71) m=0 m=0 

71). 


T These are special cases of Sonine’s second finite integral. See ref. 5, p. 376, section 
12.13, 
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and also Y = Bili—k?)! Bi —y?/D?)!, (82) 
where y = vd, Dp = wi = vf, (83) 
it follows that, for v < n, 
ny(n,v) > Be-®[I,(Y)+ BL(Y)/Y]. (84) 
For v = n, since 
(1—pyret = #21 + 48)+-0(5), (85) 
i n- 


le-B 





n| y\ n,n) 


- Be-®[y(¥)+ BI(Y)/Y],-_.—3Be-#(1 + 4B) = }Be-#(14 4B). (86) 


0 
The results for y(n, ), y(n,n—2), for example, are easily checked by direct 
calculation. 
The results in (84) and (86), written in the form 
y(n,v) ~ Be-®](Y)+ BL(VY)/Y]/n forv <n) 
1 


ke-B4 4 Be-4[1+-3B]/n for v n| 


where Yr Bi —v?/n?)!, (88) 


(87) 


may be regarded as approximate computation formulae. For example, for 
a case considered previously, n = 15, r = 5, B = 3, the following table 
shows 10® times the error in using (87)—i.e. 10° times (the correct values 
of y(n,v) minus the values calculated from the terms shown in (87)). A 
comparison of the values of the errors in Table 1 and 2 shows that for this 
case (87), up to and including terms of order 1/n, is more accurate than 
(32) up to and including terms of order 1/n? inside the curly brackets. 
TABLE 2 
i I 3 5 ; 9 II 13 15 


10% (error) 114 116 118 117 60 | —93 64 1603 


7. Passage to a continuous distribution. Connexion with the limit 
of the characteristic function 
Still pursuing a mathematical, and not a physical, investigation, we may 
now consider the limit of the discontinuous distribution as a continuous 
one, except at y | D, where the particles heap up, with a finite fraction, 
ke-®, of the total number at each of these points. Since, in the discon- 
tinuous distribution, the particles are at nd, (n—2)d, (n—4)d...., we must. 
in the continuous distribution, take 2d as the differential, dy, of the dis- 
tance y from the origin; in (84) or (87) we therefore replace 1/n by dy/(2D). 
and the proportion of the total number of particles between y and y+dy 
from the origin after a time 7 is 
B 


pee | htY)4 BAG) dy, (89) 
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Including the fraction }e~% at each end, } of the particles must be on each 
side of the origin, so we are led to the formula 


D 
Ee BY {1(Y)+ BL(Y)/Y} dy = 4(1—e-*), (90) 
ie. with y = Dsin x, Y B cos x, (91) 
( | J,( B cos x)eos y+ I,( B cos x)] dy = (e® —1)/B. (92) 


Similarly, from the formula for y? in (74) we obtain 


e® 1 sinh B 


BB BS (99) 


| I, B cos x)CcOs x7 r\ B cos x) |sin?y dy 
0 

Higher moments may be considered in the same way, but instead of con- 
sidering each separately we may take the limit of the characteristic 
moment-generating) function, which also gives an alternative method 
of finding the limiting density distribution itself. To put the formulae 
into a suitable form for the limiting operation, denote by F(y) the proba- 
bility of the algebraic distance of a particle from the origin (measured 
positively to the right) being < y, and write 


/ 


f(m) ( eiuu dF s y(n, vje’"4 — M(n, ud) = M(n,uD/n), 


} ? 
where M(n,u) is given by (16). Then in the limit f(u) > f*(u), where 


sinh( B?—u?D?)!) 


94 
(B2—u2p?)t 


f *(u) f #  cosh\ B?—u?D?)!+ B 
the convergence to the limit being uniform in any finite interval for u. 
sy general theory (6) the function F'*(y) whose characteristic function is 
f*(u) is uniquely determined, and, f *(uw) being continuous at wu = 0, F*(y) 
is the limit of F(y), at any rate at every point of continuity of F*(y). In 
lact, if F*(y) is defined as the limit of F(y), F*(y) is known from our 
previous work. It has jumps of $e-% at y = +D, and 


F*(y) 0 tor y D 
_ Be-B ff... ee ee — 
Le-B 4 ™ [1(¥)+BL(Y)/Y]dy, for —D<y<D, 
2 J, 
D 
a . — ‘a 
i 3D [1(Y)+BL(¥)/Y|dy, for y >D, (95) 


. 


D 
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where Y is defined in (82). To prove the connexion between F'*(y) and f *(x) 
it is sufficient to prove either that 


f*(u) | eiuy dF* — eB cosuD + aed [J)(V)+ BL(Y)/Y |cos uy dy, 
y=—« 0 (96) 
i.e. that 
[ | I,(B cos x)cos y+ J,(B cos x) ]cos(uD sin x) dx 
0 


cosh( B?—u?D?)! | sinh( B?—u?D?)! cosuD (97 
= —_———_—_— a i _ — F { i) 


B (B?—w?2D?)! B 
or that 


27 cosh(B?—u®D*)' | sinh(B?—u®D*)! _ cosuD] 4, 
| | + ee —B 
oe plo) + BAY) 7] Seoae<d 


I é 7 r , 
“ spol} )+ BL(Y)/Y] aty=D,Y=0 


-Q0 fory > D, (98) 


each of which is the Fourier transform of the other. The first of these may 
be deduced from Sonine’s second finite integral (ref. 5, p. 376, section 
12.13), and the second from Sonine’s discontinuous integral (ref. 5, p. 415, 
section 13.47). 

This section may be closed by a final remark about the continuous 
distribution (89). We have taken the limit when » — 0, 7 + 0, with T 
(=r) fixed. As 7 varies for the continuous distribution it will be the 
constant A (= limz/(1—c)) that will be independent of 7. In the limit 
c>1, and B = 7/2A. If we make this substitution, the formula for y’ 
in (74) becomes (2) or (22), and for large 7’, 


y? ~ 2Av*T. (99) 
It is easily proved that 
spo? |a+ a 1 ~ (2ry*)-texp[ —y?/(2y?)], (100) 


with the above value for y?, so the distribution tends to the Gaussian; and 
it is easy to verify directly that the characteristic function for the variable 
y/(y?)* tends to exp(— 4u?). 
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8. The governing differential equation as the limit of the difference 

equation. Solution of the telegraph equation, without leakage, 

for an instantaneous plane source 

Let us now consider the matter in a rather more physical fashion. (To 
begin with, change the notation slightly, and write ¢ instead of 7' for the 
total time of the motion.) If there are 2.N particles altogether, we see that 
they spread out symmetrically from the origin, NV to the right and N to 
the left, with a ‘wave-velocity’ v, and that those to the right heap up at 
y = vt, where there are Ne~?4, that there are none at y > vt, and that in 
) << y < vt they are spread out with a line density No(t, y), where 
etl2AT , t TY) 


a(t, y) I(¥)+ (101) 
2PAv ‘ 


(vi? y?)! 


and r oy 


(102) 
We may further consider the process involved by deriving its governing 
differential equation as the limit of the difference equation (12) for y(n, v). 
The differential of time, dt, will correspond to the interval 7, and, as 
previously explained, the differential of distance, dy, to 2d or 2vr = 2v dt; 
ewill be replaced by 1—dt/A and p ( 4(1+c)) by 1—dt/2A, so, apart 
from the discontinuity at y = vt, o will satisfy the limit of the equation 


o(t+dt, y) 


(1 sa )ley kdy)-+-o(t, y+ 4dy)|— (1 — 4 Jetta). (103) 


If we neglect terms of a higher order than (dt)?, this becomes 


e 
Co l o*o 








o+ dt + (dt)? 
ot 2 ot? 
dt l Co dt Co l eo 
] 20+ (dy)? ( . ees, dt)?|, (104 
| 9A} | "4 aye? Aa +3 oe | O%) 
2 | sf c2 
ie. (: r+ (dey = May), (105) 
ot? A at oy” 
so the limiting differential equation is 
Co | l do poe. (106) 
et? ' A at oy? 


This is the telegraph equation, without leakage, and with v? = (LC)-, 
A= L/R, where L, C, and R are the self-inductance, capacitance, and 
lesistance per unit length. The equation is satisfied by the potential, and 
also by the charge per unit length, of a cable; the potential and charge 
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per unit length are proportional to one another in the theory. It is easy 
to verify directly that (101) is a solution of (106). It is clear from the 
process by which it was obtained, and it may be verified a posteriori, that 
the problem solved is that of the ‘diffusion’, according to the telegraph 
equation, of an instantaneous source at the origin; since the ‘diffusion’ js 
in one dimension, the source may be called, by analogy with the corre- 
sponding problem in the theory of the conduction of heat, an instantaneous 
plane source; if we have a straight cable, infinite in both directions, and at 
time ¢ = 0 introduce a charge of 2 units at the section y = 0, then at any 
time ¢t there is a charge e~?4 at y = vt and at y = —vt, together with a 
line density of charge o given by (101) between these two sections, and 
no charge for |y| > vt; A is L/R and v is (LC)- as above. 
This solution is analogous to the well-known solution 


l ( y? bi 
exp! —- (107) 
(arKt)? At 


for the heat content per unit length in an infinite medium of thermometric 
conductivity «, when an instantaneous plane source of strength 2 is intro- 
duced at y = 0 at time t = 0. (The equation of heat conduction is known, 
and easily shown, to be derivable by a limiting process, similar to that used 
above, from the differenee equation (12) when the motions of a particle in 
two successive paths are uncorrelated, i.e. c = 0, p = q = 4; the thermo- 
metric conductivity « is }v27 = 4rd; as n > 00, tr > 0 so that nz = t is 
finite, but v must be allowed to tend to infinity with n! (in fact v? = 2«n/t) 
and d ( vt) > 0 like n-!, so that nd becomes infinite; the differential of 
the time dt, in the limiting process, is of the order of the square of the 
differential, dy, of the distance y. The limiting process in the case con- 
sidered here requires v to remain finite, and to be fixed, instead of tend- 
ing to infinity; nd remains finite, and dt and dy are the same order; the 
root mean-square distance from the origin in the discontinuous process is. 
in our case, proportional to n for large n, instead of to n! as for ¢ = 0. 
For more general discussions of the passage to a continuous system from 
the problem of the random walk, reference may be made to refs. 7 and 8. 
The case considered here appears not to have been previously discussed.) 


9. A solution of the telegraph equation in various forms, and 
mathematical deductions 
A more usual problem in the theory of the telegraph equation is that in 
which a semi-infinite cable (i.e. for practical purposes, a cable of length 
greater than vt for the greatest time ¢ considered) has one end at y = ", 
and initially the potential V, the charge, and the current are all zero. At 
time t = 0, V is raised to unity at y = 0, and it is maintained at unity at 
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that section. It is required to find the values of V along the cable at 
subsequent times. The analogous problem in heat conduction is that of 
finding the temperature # in a semi-infinite rod, when # = 0 initially, and 
$is raised to unity at y = 0 at time t = 0, and is then maintained at that 
value at y 0. The solution of that problem is 
y y|{2(xt)*} 

> PT 


l y* 2 : 
i=] exp| J dy = e-* dé i—erf —4 
(art)? 4 ict \7 2(«t)? 


. 
0 0 


(108) 
and the required solution of the telegraph equation is related in the same 
way to (101) as (108) is related to (107), account being also taken of the 


finite concentration, e~!?4, at y = vt. The solution is 
' 1,() | 
om FT t 1,(Y) . 
, 1 aa a 
J l a(t, y) dy l Av | [,(Y¥) + vA - dy fory < vt | 
0 0 


| 
0 for Yy > vt, 


(109) 
where Y is given by (102). Since éc/éy = 0 at y = 0, V satisfies the same 
equation (106) as o. Also, as we have already substantially pointed out 


in equation (90), we must have 


| o(t,y) dy = 1—et4, (110) 
This equation reduces to (92) with ¢/2A in place of B; its verification is 
therefore included in the verification of (97) (with w = 0 for this particular 
case); it is also verified quite simply by expanding the integrand in powers 
of (t/2A) and integrating term by term; cf. ref. 5, section 12.11, pp. 373, 
374.] Hence lim V = e-2A (111) 
/ vt —0 
and V has a discontinuity, with a jump e~?4, at y = vt. 
If p is used now to denote a Heaviside operator (its previous use as a 
probability is not required in this section), the operational representation 


of V is ; 
Yy 2 pP 2 

24 . 112 

‘(p ‘) | (nt?) 


and for y vt Jeffreys (9) gives a solution in series 


a a W mn r 
7 e-t[2A VY)+2 > Y)}. 113 
J fi 4 (>) [,( ] ( ) 


n=1 


where w.%9 *_(@—" (114) 
2Av Y vt-+y 


V = exp 

















150 8S. GOLDSTEIN 


Although it leads us away from the original purpose of this paper, it is 
of some interest to discuss mathematically the connexions between the 
different expressions of the solution. In the first place, from (112), 


c+ia 
, 1 f¢ y[{.. , 2\4) dz 
, = — exp/et—¥ 224 | (c > 0). (115) 
2m | v A | z 
c—iz 
V = Ofor y > vt is an immediate deduction. For y < vt, the contour may, 


by Jordan’s lemma, be transformed into —oo, (—A+,0-+-), and may be 
taken to be the lower and upper sides of the negative real axis, indented 
by small semicircles at —A, and completed by a small circle at the origin, 
the radii of the semicircles and the circle being finally made to tend to zero. 
In this way, and with z —ax = —(l+ cos €)/(2A) on the real axis, (115) 


becomes Pa 
24 
- é 
V=1— | exp 


7 


tcos € 


2A 


dé, (116) 





ysin€]  sin€é 
2Av 








]+-cos € 
0 


and the identity of (116) with (109) is established by integrating (79) and 
(80) with respect to y from 0 to y (with B = t/(2A), AB = y/(2Av)). If 
in (115) we make the substitution 


us 





hice tenth oe (117) 
24 2AW 4¢ 
in the integrand, the integral becomes 
14 7. yLIiWw 
"ae ee | walt. C aW de (118) 
2a | 4¢/C—43W ¢ 


and, for.y < vt, we may take the contour of integration as 
-0, (—4W+,0+,4W-+). 


We now introduce the Lommel function of two variables U, (ref. 5, 
sections 16.5-16.59, pp. 537-50). From ref. 5, section 16.58, p. 548, 


1wW+,0+.4W 4 


(—34 2 ) 
ne eee y) (w/t de 
U, (iW, iY) = mpl 4) es =: 
n(tW,cY) | exp a) 1—i We £ a 


D 
a7 
© 


(0+) 
. 1 re 
Since also L(Y) = = F | exp(+ i ° (126) 
“7 J 


(a Bessel-Schlafli integral), it follows from (118) that, for y < vt. 


V = e24(2U, (iW, iY) —2iU,(iW, iY) —h(¥)]. 
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Jeffreys’s series (113) follows at once from 


a "\ 2m ne "\ 2m+1 
UiW,iY) = > ef hy(Y).  UiW.AY) =i > (F) hsh¥ i 


m=0 m=0 (122) 
= 


ref. 5, section 165, p. 537). 
That (109) and (113) are equivalent may also be simply seen by differen- 


tiating the right-hand side of (113) term by term with respect to y, and 


using the recurrence formulae for the J, to show that the derivative is 
[1)(Y)+4(W/Y+Y/W)L(7)}- 


The following formula follows from the — of (109) and (113). 


[,(zsin@)+2 > tan? 367, (zsin@) = 1+2z i [ /,(z sin x)sin y+ 4,(z sin x)] dx. 
n 1 
(123) 
[Since e? = I(z)+2 > I,(z) 
n=1 
(92) is the special case of this when 6 har. | 
[It is also of some small mathematical interest to remark that 
> g(vt—y\*” 
na(7 9)" 17) (124) 
vt+-y 
is the solution for y < vt when the end condition is V = e-?4J (t/2A), 


whether » be integral or not. The operational representation of this 


| p y l wre: yl. , p\t 7 
2 7» 4 | et e >. 4 ~ > he me 12: 
(2A) e va) [ 2A (p I | — ‘(» A me 


that (124) is the interpretation of (125) may be established by proving 
either that 


solution is 





l 
tan” 46 [ (x sin @) | exp(x cos €)cos(x cos @sin €+-n€) dé— 


exp(—a cosh +2 cos @sinh E—n€) dé, (126) 
or that 


[ exp| — Z(cosh €— 1) |Z, (z sinh €)tanh" $€ sinh € dé 
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of which one is deducible from the Mellin transformation of the other: 
(126) may be deduced from the formula (5) on page 177, section 6.21, of 
ref. 5 (due to Sonine) for re[a(1—cos @)| > 0; (127) may be deduced inde- 
pendently, for re(Z-L-z) > 0, by integrating term by term the expansion 


1,|2(u?+-2u)*] 7” S zm—n I wi (128) 
1 


2?(u?+ 2u)in m!u" 
( 


after multiplication by e-“2u”, and using the formula 
« io) 


4 a—(a?—b?)|" 
e-@"T (bu) du _| ( é‘ ’) | (129) 
L b"(a?—b?)! 
0 
(ref. 5, p. 386); the expansion (128) is a slightly altered form of Lommel’s 
expansion (ref. 5, section 5.22, p. 140). 


10. Further developments. The problem with leakage included 

Certain generalizations of the preceding work are suggested. One of the 
most immediate and obvious, which is not uninteresting, is to generalize 
the work to include leakage in the final form of the telegraph equation. 
The whole discussion has started, of course, as a generalization of the 
classical problem of ‘the drunkard’s walk’, and the further generalization 
required is this: to find the probability, after he has had time for a fixed 
number of steps, that a drunken man will be found a certain number of 
steps forward from his initial position if he manages to keep moving in 
a straight line (if he does not disappear) and there is a probability p at 
each step of forward motion and a probability g of backward motion, but 
also a non-zero probability, 1—p—gq, after each step that he will, in fact. 
fall down a man-hole and apparently disappear altogether. In terms of 
the motion of a large number of particles, with which we started, we must 
dispense with the condition p+q = 1, and allow a non-zero probability 
1—p—gq that a particle will escape from the system at the end of a 
partial path. We may define a(n, v), B(n,v), y(m,v) as before, but we must 
now add that they are to be taken at time nr— 0 (i.e. before the disappear- 
ance of a certain fraction of the total number of particles at time nz). We 
still write c = p—q, so that c is the correlation coefficient between the 
directions in two successive partial paths for those particles that do not 
escape from the system; we also write r = p/q as before. Then many of 
our equations are unaltered—e.g. (4)-(10) inclusive. Equation (11) is 
replaced by 


x(n, v-t 1)+-B(n, v—1) (p+q)y(n—1, v), (130) 
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ind (12) and (15) by 





y(n, v—1)+ y(n, v+1)|—(p?—¢?)y(n— 1, v), (131) 
ind M(n+1, u) 2p cos uM (n, u)—(p?—q?)M(n—1, u), (132) 


y(n+ 1, v) Pp 


respectively ; (14) is slightly altered to 
M(1) = cosu, M(2) = q+pcos2u = p+q—2psin*u, 
ind the solution for M(n,u) is still of the form (16), with a, and a, still 
given by (18), but with 
k, l v me = te l . = ss i} 
2(p+q) (1—r* sin*u)? 2(p+q) (1—r*sin*u)? 
(133) 
If we express the previous values of M(n,u) and y(n,v) in terms of 
p/q. in addition to n and u, and to n and v, respectively, then the new 


values are the same as the old values multiplied by (p+gq)"-!. We now 
have, for example, ‘ 

yd y(n, v) (p+q)""}, (134) 

n e r 1 n\o ee 

and > v*y(n,v) (p+q)"*|nr—4(r? yf) ( ) \| (135) 

t nm T+ 
corresponding with (20). 
If we put | — = A, (136) 
oa 
is before, and also = G, (137) 
1—(p+q) 
ind let n + co, r > 0, with nr = t, as before, then in the limit the mean- 


square distance from the origin is given by 


2t 2 sd 
ye = vee | —{1—eWld-i@], (138) 
] A 1/G (1 A ] G)* 
corresponding with (22). 1/G@ is the proportion of the charge that leaks 
way per unit time. 
If we now form a differential equation from (131) in the same way as 
106) was formed from (12), the result is 
eo l l\ do o > a ‘ 
} eee te, (139) 
1G) et AG oy” 


ot* 
which is the telegraph equation with leakage, with v? (LC)-! and 
A = L/Ras before, and with G = C/S, where 1/S is the leakage resistance 
per unit length. Since the charge per unit length is CV, where V is the 
potential, and the leakage current per unit length is SV, this is the result 
we should expect for G. 
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Equation (139) may be written 


a’ ( 1\ do’ 2 00’ e-llGg’ 
+{—— = Vv —_, = 2-7-8, 
2 A G) @t oy? lead 


ot 
and we may easily convince ourselves that not only (138) and (139), but 


a9 


C 


all the limiting results, are obtained from the previous results by replacing 
A-! by A-!—G~-! and multiplying by e~/%. Thus (101) is replaced by 


1/1 1 t{/1 1 . tf1 1\ L(Y 
ay) = sel 7 a)exr| —5( a)| {tor +54 7m al aa me) 


where Y is now given by 
‘ 1/1 1 \? y 
c* = ——| |e — 142 
4\A G Yr i 


in place of (102). (It should be added that, in order that q should be posi- 
tive, 1/A—1/G must be positive.) The finite concentration at y = +t is 
now exp| —43¢(1/A+-1/@)]. 

The result corresponding to (109) is the solution of the problem when 
V =e at y = 0fort > 0, and is obtained by replacing A-! by A1—G" 


"4 as before. To obtain the solution of (139) with 


and multiplying by ¢ 
the end condition V = 1 at y= 0 for ¢ > 0, we must start from the 
solution of (106) with the end condition V = e’@ at y = 0 for t > 0, and 
then replace A-! by A-!—G-! and multiply by e-/%, as before. The 
required solution of (106) may be found from that given in (109); if Vi(t,y) 
denotes that solution for y < vt, the solution of (106) with the end con- 
dition V = f(t) at y = 0 fort > 0 is 
t 
V =foyhity) + | £e—OK(é.y) dé. (143) 
ylv 
where f(t) is df/dt. It may be added that, with leakage present, the solu- 
tion for V = 1 at y= 0 for t > 0 may also be expressed in terms of 
Lommel’s function U,, of two variables: it is (for y < vt) 


V exp —5( 4 a) }(Caem.«Y) + U,(iW,, i¥)—iU, (GW, iY) 
~iU,(iW,,i¥)—L(¥)], (144) 


.. pw fl. 1\V y — if. ly iy , 
W here W, = 5) (4i+ zi) (4), W, =z 5) (a 7) ( ) (145) 


This solution holds whether A-1—G-1 is positive or negative; Y must be 
taken as the positive square root of (142). Jeffreys’s series for this case 
is immediately deducible, but I have not pursued the matter further. 
Other generalizations lead to higher order difference and differential 
equations. I have, for example, briefly considered the case when, in 
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addition to non-zero probabilities that a particle will move forward or 
backward with a velocity v, there is also a non-zero probability that it 


sa | will come to rest (the swaying drunkard), together with a non-zero proba- 
| bility that it will start again from rest with velocity v. This leads to 
ut ‘ . . » . . . i . . . . 
, third-order hyperbolic differential equation in the limit, containing 
cing ; “ > : ‘ 9 9 [0.46 ° : 
@o et, C2a/et?, Ga/et, Ba/ey*et, C2a/ey*, the highest terms being 
y { 
a/ ot? v? Bo oy et. 
14 , : . , —" 
| (ther problems arise by introducing non-zero partial correlations between 
the motions in non-consecutive intervals. 
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r. P(a,b;c;2) is the hypergeometric function, defined by 
rential | . a.l u(a+1)b(b+1 
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(28) may be written 


rn—3(7— 1) n—v ] 
y(n,v) = F(— g tty +451; 5)4 


(r+1)" 1 

n—3 typ 

‘. a (F(—" ry. Bt a ‘] 
y r+ n & z 2 


for v < n. It is, in fact, a known result that the solution of a difference 
equation such as (12) may be expressed in terms of hypergeometric func- 
tions. In accordance with the notation of section 6, we now put r = n/B, 


— 


v = Kn, and let n > 0. From a result given in ref. 5, section 5.7, p. 154, 
it is known that 


9 


] rs 
li FiA,u;m+1; 
im — ( pesm i 


Amo Mm! 
ll 


= 1,(¥)/(4Y)", 
fA 


where Y = limz. The limiting form of y(n,v) for v < n, given in (84), 


follows. (The limit given in (86) follows at once from (27), since 


p= r/(r+)) (1+ B/n)-1.) 
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SUMMARY 


Rayleigh’s problem, in which a plane is started to move parallel to itself through 


scous fluid, is considered for the case when the fluid is compressible. An approxi- 
ite expression for the pressure-time variation at the plane is obtained and the 
nitial motion is considered. It is shown that the initial pressure variations deter- 


ned, as here, from the classical equations are so rapid as to require reconsideration 
the light of molecular phenomena, though there is little doubt that the entire 


ressure changes are limited to an interval, short by ordinary standards, and 


f 


governed essentially by the time relaxation. This conclusion in itself is thought to 


ye of some interest in the theory of viscous compressible flow, and its implications in 


relation to the steady flow along a semi-infinite plane are considered. 


1. Introduction 

THE motion set up when an infinite plane, immersed in viscous incompres- 

sible fluid, is set moving instantaneously from rest in a direction parallel 

to itself (and thereafter maintained in uniform motion) is one of the stan- 

lard problems of incompressible flow. The motion is everywhere parallel 

to the imposed motion of the plane, and the velocity is given by 
U'(1—erf y/2(vt)*), 

where U’ is the velocity of the plane, ¢ is the time, y is the perpendicular 

listance from the plane, and v is the kinematic viscosity. 

In the present paper we shall consider the corresponding problem when 
the fluid is compressible, and we shall, in fact, limit attention to a perfect 
gas. Compressibility, of course, considerably complicates the problem 
since the temperature variations, which occur through dissipation, induce 
density variations which in turn give rise to a component of flow at right 
ingles to the plane. We shall limit attention to two aspects of the problem: 


an approximate solution for the pressure—-time variation at the plane 
when the Mach number is small and the Prandtl number o 3; 
li) a representation of the initial motion for the full equations when 
temperature variations of the coefficient of viscosity are neglected. 
The approximate solution referred to in (i) is obtained by linearizing the 
equations in all save the dissipation terms, but is an improvement on the 
analysis of Lagerstrom, Cole, and Trilling (1) (hereafter referred to as 
L.C.T.), who neglected conductivity and did no more than obtain a re- 
peated integral in complex form for a general value of the dissipation. 


({Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 2 (1951)] 
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We shall allow for conductivity in the special case when the Prandtl number h 
is } and shall derive, at the plane, the particular solution associated with ae 
the actual dissipation as a single finite real integral, which will be tabulated, Fir 
It should be said at the outset that, for air under standard conditions. 
the solution derived below on the basis of the classical equations gives 
significant pressure variations in times of the order of 5 x 10-1 sec. Modi- 
fications of our results in the early stages must therefore be expected in 
the light of molecular considerations. Further discussion of this feature 
of the solution is reserved until section 5. and t 
2. Equations of motion 
With Ox measured parallel to the direction of motion and Oy perpen- 
dicular to the plate the equations of flow are in the usual notation* 
Du o[ eu 4 
= AL 9 
P Di oy ‘ oy 
Dv op 4 ef[ oa 5 
Dt ey 3 ey\" ey)’ “) 
op @ 
4 —(pv) = 0, (3) 
ot oy 
D é é 
where = —+v_—. (4) | where 
Dt et! ey \ : 3 
iia ; and J 
lo these must be added the equation of state 4 
An 
p = RpT (5) 
and the energy equation 
DT Dp é [,,eT éu\? 4 /av\? 
pJc, - : = J B k 7 +HII> + 7 (9) |, 
Dt Dt ey\ ey oy 3 \ey The b 
From (3), a function y% exists such that and ti 
Cn yp Py= 2 (7) 
Ps Y Ps ot 
where here and elsewhere the suffix s is used to denote some standard con- 
dition of the gas which will, in fact, be the undisturbed state. 
We may, then, in von Mises’s fashion, transform to x and y as indepen 
dent variables as indeed Illingworth (2) has done-in the boundary layer | 3. py 
approximation to the present problem. Then we have If J 
(=) - (5) _ 2 (5) (8) | small 
ot}, ot}, ps \ep), if we» 
(=) _P (; ) (9) p 
Cy}, — ps \O], 
* All x-derivatives are clearly zero. the 8q 
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- D C 
poe so that : ( : (10) 
with Dt ly 
ted, Finally putting 
ions 
ives us V8’, t ar. w= Uw v= a, 
; a ae 
lodi ‘ ; 
din p PsP’; p PsP's T 1,7" ee) 
ture 
ind then dropping the dashes we have 
Cu C ou 
. (12 
ot al ps) 
‘pen- | 
| ~ 4 Das l y 
ofa) MS (3) 
ot 3 ob ous Y ous 
‘ (? ov (14) 
| ct\p Cus 
2) | 
esol (15) 
| 
7 hs 2 4/dap\2 l < oT : 
| tito a etl) 
ct p ot ous 3 Ors o Oy op 
4) | where o has been written for uc,/k, y is the ratio of the specific heats, 


| ind 
| 


M U /a,. 


An alternative form of (16) obtained by using (14) is 


ou\? w\2] lal af 
ee mu Ph. a) ri J \ a7 
ty) 0 , (5) 3 4 o ap \Y? ap ati 


lhe boundary conditions are, assuming the plate to be thermally insulating, 


ov 


1 oT y—1 
, at | 


Cc us 


nd taking conditions in the gas at rest as standard, 





i | int ent Sn oe ~ = 0, 
| Oy 
dcon- | 18 
da con ” 0. 99 0. p l. p 1. ( ) 
ch latt 0 and as f > oo. 

depen 

reT e ° : 
 layeT | 3. The linearized solution for small M2? when o = 3 


4 
| If M? is small, the variations of pressure, density, and temperature are 
s) | mall and so is the velocity component normal to the plate. Therefore, 


+ 


i we write 


; p=1+p’", p=1+p’, T=14T", p=l4+p’, v= 140, (19) 


| ‘he squares and products, of the dashed quantities may be neglected for a 
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first approximation and we have from (13)-—(16), again dropping the dashes. 


ov 40% lép 


7 — 9 <« a ‘ ? (20) 
ot 3 Op? sy Ob 
ec 0, (21) 
ot ° ods 
nie Ys ses l = ea. ] = ‘ 
ot y a o ys” 
p= p+T, (23) 


where © is the dissipation function which, to a first approximation, is given 
by 


® = (y nar(=). (24) 


The velocity wu itself is O(1), but since it affects (20)—(22) only through © 
which contains the factor M? it is sufficient to determine the leading term 
only in wu, and this is clearly given by 


‘ a9 
ou oru 


~y — naDe? (25) 
ot Cys? 
so that u = 1—erf(s/2#*), (26) 
| _1)M , 
and therefore Oo — ly ; eet (27) 
to the order required. 
The boundary conditions on v, p, p, T' are now 
; oT 
(i) v= 0, =0 wheny = 0, 
us Y (28) 
(ii) v,p,T,p>0 asp>oo. 
Now it follows from (21) and (23) that 
oT : op, wv (29) 
et at | ep’ 
and from (22) that 
(; _2 ier y—8 ep oD (30) 
6t a Oy?) at y ob at 


Therefore, inserting the value of @7'/ét from (29) in (30) it follows that 
ra) i e Op a. e ov _ y—1 ep | 6D (31) 
to ob) ot | \et a ep?) ep y ot | at 

Comparing (20) and (31), it is evident that ¢ = ? is a specially simple case 

and, since this value is so close to the observed value 0-7 for air, we shall 
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take advantage of the simplification offered and limit attention to this 
value. 


Then substituting from (20) in (31) we have 


é 40°7\0p 10p y—1ap oa® 
» 5) : 2 - ‘ = ‘ ? (32) 
\ct Seb) et yap y oa ot 
4 O'p Cp =p om 
that is, —t 4 #. A —- p —y—. (33) 
3° Ow? ct Onb* ot ct 


This equation is exactly of the form obtained by L.C.T. (1) who neglected 
onductivity altogether. L.C.T. (1) considered the solution of (33) for a 
general value of the right-hand side only and used influence functions to 
express the solution as repeated integrals. We, on the other hand, shall 
make use of the known form of ® and shall consider the solution of that 
specific equation. 


The boundary conditions on p may be determined as follows. Firstly, 


since v 0 when us 0, it follows from (20) that, at 4 = 0, 
4 = ] op (34) 
3 ous? y Ous 
Also, from (29) by differentiation 
aT 2 52), 
at ETc (35) 


Otords ots ' Cys? 


Now, since we are considering a thermally insulating boundary, é7'/éys = 0 
0 and, since this is true for all time, é*7'/etés = 0. Therefore, 


when & 


from (34) and (35), an 8 op 
. 
. 4 P_9o (36) 
Ctous 4y Cus 
when us 0 for all a 


That is, ( a Ae-*i4y, (37) 


But (ép/eb) -0 when t + 0;* (@p/és),-9 must therefore vanish for all 


time. Hence our boundary conditions on p are 


cp 0 
ot 


when f 0 (4b > 0) and when & > -, 


0 when wd 0 (t 0). 
Hence, with a convenient change of variables, 
T 3t/4y, x 3s 4y, (38) 


See, for example, the initial motion in section 4. 


M 
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we require to solve the equation 


23 a2 22 © 
Gp _ep ¢ p Bp C I 4yx?/37)_ (39) 
Ox*0r Ox? «= Or® 2 Or\r 
where B = —2y(y—1)M?/z, 
with the boundary conditions 
cp 
p= Z = 0 whenr 0 (xy > 0) andasy>« 
CT 
‘ (40) 
cp 
a = 0 when y= 0(7>0). 
ox 


This is a problem admirably suited to the application of Heaviside’s 
method since the right-hand side can be expressed operationally in terms 
of Bessel functions.* To avoid confusion we shall henceforward write P 
for the excess pressure and p for Heaviside’s operator @/ér. Then, since 
67 ’ l cee 
pK p*x h') exe (41) 


) = 
aT 


(where K, is Bessel’s function of the second kind of order zero) we may 
put the right-hand side of (39) in the form 

Bp? Kl x(8yp 3)! |. (42) 
Hence the operational form of (39) is 


a9 > 
o* , oI ‘ 
—.m*P = Bm?K,| pay], (43) 
C : 
where m2 = pi(p+1) and a = (8y/3)!. (44) 
Variation of parameters leads to the operational solution 
: : 
P —Bm|coshmy | e~™K,(ptat) df+-e-"x | cosh m£Ky( p'at) dl}, (45) 
" xX 0 a 
and, in particular, the excess pressure P,- at the wall is 


4 


w -Bm [ e-™ Ky( pial) dl (46) 
0 
—Bn [ e-"#K,(u) du, (47) 
0 
where n = m|(ap)? = [ p/a(1-+p)]}. (48) 


Hence (see Watson’s Bessel Functions, p. 388), 
4 = —(Bneos—!n)/(1—n?)!. (49) 


* For an account of this method and its relation to the Laplace transform the reader is 
referred to Jeffreys and Jeffreys (3). 
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The expressions (45) and (49) have then to be interpreted by means of 
3romwich’s integral, so that we have, from (45), 


t 











8 , e’t } A q a y 4 
P ‘cosh | —“* eA +" K(AaL) dl4 
2ri J (1-+A)!] (1-+A)#] . sell 
x 
e—AxiL +A)" | cosh . K (Mat) dt dr. (50) 
. (1-+-A)! } 
0 
As L.C.T. (1) have shown, the exponential essential singularity at A | 


| makes further progress difficult except by way of asymptotic expressions 
| . . ay . . 
| which are clearly insufficient when, as here, the result has to be integrated 
| again. As they also remark, the functions represented by integrals of the 
type 
| ‘ ; 1 b(A) 
r a d 


eAt—Ax\(1 


271 A * (51) 


where (A) is reasonably well behaved, will have to be studied further 
before the determination of conditions throughout the field of flow can be 
ompleted. 
| Fortunately the form given in (49) when y = 0, similarly interpreted, 
| has only a logarithmic singularity at A | and can be dealt with satis- 


factorily. For, in virtue of (48) and (49), we can write 


3p v2 + (a?—1)p]4 . 
Py i 5) tan | ai (52) 
| (a"—1)p- xv? |! p 
ind by Bromwich’s integral the interpretation is 
3 er y? +-(a2—1)A]4 a 
Ey —, tan- V2 ay. (58) 
27i J AA (a®?—1)A+a2?| A 
The integrand has branch points at A = 0 andA x? /(a?—1) —b say, 
ind A |. The last named is in fact a logarithmic singularity; for 
tan-!z -log a 
21 iz 
} and has branch points (like logz at z 0) when z +4 (i.e. 22 = —1) 
which corresponds in our case to A i. 


The usual arguments suffice to show that the contour of integration can be 
modified to the two sides of the negative axis indented at the three branch 
points. The integrals about each of the parts of the three indentations tend 
to zero with the radius of the indentation. Hence we are left with the 


integrals from the various straight portions. 
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The three branch points make the evaluation tedious without being in 
any way novel and the analysis will not be reproduced here. The results 
obtained are as follows: 

(1) The contribution to the integral from the two sides between —b and 

2% and between —1 and —b are found to be of the same magnitude and 
opposite in sign and therefore cancel. 

(2) The contributions between 0 and | give 

F 1 sin-'b-! 

: = ‘ 
a a en | ebrsin’ d0. (54) 
2(a2—1)* J pi(b pe)? (a*—1)? . 


0 0 


) 
UL 


The pressure therefore rises instantaneously at time zero by an amount 


- . se Sy—3\} 
(Py h=0 Vly 1) ; 3) sin : . "9 (55) 
7 Sy 3 8y 


which may be ascribed to the infinite initial dissipation. This result is 
confirmed in detail by the analysis of the next section where the manner 
in which the boundary conditions are satisfied is made apparent. 

A few values of the excess pressure Py, (which is, it should be remem- 
bered, effectively expressed as the ratio to the undisturbed pressure p,) 
are given in Table I, for y = 1-4, as a function of a2t/v, (where ¢ is the 
actual time). 


TABLE |] 


1/5 P.M? 

0°0000 o°22I 

1*°306 o'17i 

°733 07138 

5°407 or1ol 

10°93 0:070 

21°87 0°049 

F v, \} 
Asymptotically iy o22sar{ 7) (56) 

as 


4. The initial motion for the complete equations, temperature 
variations of ». being neglected 
We return now to the full equations of motion set out in (12)—(16). The 
method to be developed here would, in fact, be possible for any given 
temperature variation of u, but the details are considerably simpler when 
pw is constant and we shall restrict our discussion to that case. 


The solution of the corresponding problem for incompressible flow sug- 


gests the use of a new variable 





yb t}. (57) 
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in and we shall for convenience associate with it as the second independent 
ts variable , _ 
7) t?. (58) 
: . = C l ( a C Lé ‘= 
nd =~ Then -| 7 £—| and —: —— (59) 
ot 2n*\ On ag ous n O§ 


nd 
| Itis then possible to satisfy the equations by expressing u, v, p, and 7’ as 
power series in » with coefficients which are functions of € only. All the 
series save v are in even powers of €; v is in odd powers. The details of the 


algebra, which is tedious, are omitted. Suffice it to write down the equa- 








4 : — = : ‘ - , 
tions satisfied by wp», Po, Po, 7) the terms in the various series independent 
of , and by v, the coefficient of 7 in the series for v; they are 
in | Po ga (6C) 
5 Po 1, (61) 
Ug t+ téuy = 0, (62) 
nel is rh le oly 1)M*u,?, (63) 
» 
Q <,! _ = , > 
m U1 $l\— Py ~ Po (64) 
P d 
t om“ ° ° ° 
the where dashes denote differentiations with respect to &. 
Had we not assumed pu constant we should have obtained (60) and (61), 
but (63) and (64) would have been replaced by simultaneous equations for 
u, and p, rather than, as here, equations admitting of separate solutions. 
' It is not difficult to write down the equations for the later coefficients, 
but their solution has not been attempted and so they are not reproduced 
here. There seems no reason to doubt that a solution on these lines could 
be developed, but in the light of the subsequent discussion it scarcely seems 
called for. 
The solutions of our equations are 
18) 
uy = 1—erf 3é, (65) 
c 
ure (4 | . y 2 : , 
py) = T, = 140 M2 | e-ottty | eo-2yn*Ay dnd 
The ’ 
Ve ' | 1 . - . 2y a Ss ae 
her ( 1)M?) ‘ | e-95°/4Y erf *tidl (66) 
lar(2 oc) 4 4y 
sug P iM w - 16 P 
Yes « ¢ e-ts >" a - 
v; é é | | | eo—2y)n?* 4 dy . | eon 16 dy di. 
70 tc) = 4 


(67) 
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We may notice that, at the wall the integral for p, can be evaluated to give 


2 at oe 2y—o\} 
vy 1) }*tan— : f 
7 2y 0; Oo 


which, as already anticipated, agrees with (55) when o 


(Po) n ] (68) 
* and provides 
a useful check. 


With o 


1-4, the values for air, we have 


LQ 2155M2, 


0-7 and y 


(Po)n I 


Table 2 gives the variations of py and v, with &€. 


TABLE 2 


(Pg— 1)/.4 ay Ve 

O'215 00 
i775 ‘O214 
I O*1032 "0355 
| 455 t 215 
9163 005 
} O 49 9025 
5°55 ¢ O1 "0004 


It is, of course, not possible to determine precisely the range of validity 
of the solution obtained without determination of higher terms in the 
series. All we can say is that it will be valid for sufficiently small times 
and that these times will be measured in relation to v,/a? or v,/U*. It may, 
in fact, be shown that it is always the smaller of these two quantities which 
is the relevant comparison time. 


5. Discussion 

Let us consider, in the first place, how the results we have obtained fit 
into the pattern of any complete solution of our equations. It is very likely 
that, as Carrier and Lin (4) found in their considerations of the incompres- 
sible flow past a semi-infinite plane, the flow can be split up, effectively, 
into four regions, overlapping in parts, as follows: 

(i) a solution useful at small times; 

(ii) a boundary-layer solution valid for large times and y = O(t'); 

(ili) a region of effectively inviscid flow; 

(iv) an intermediate region in which conditions are complicated and in 

which numerical methods of solution seem to be the only possibility. 

(i) is the solution given in section 4 above; (ii) has been given by LIlling- 
worth (2). (iii) and (iv) have not been attempted here, though considera- 
tions similar to those of our section 3 could possibly be used to extend (iii) 
and join it on to the other flow regions. 
It is interesting to notice that, in the initial solution, the density remains 
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unaltered and the pressure variations, in consequence, are identical with 
the temperature variations. In contrast, in the boundary-layer solution 
the pressure is effectively the undisturbed pressure, whilst the density 
is then the inverse of the temperature. The velocity normal to the plane in 
the initial solution is proportional to # and rises from zero at the plane to 
, maximum and then falls off to zero. The density variations determined 
in the boundary-layer solution imply a velocity normal to the plane propor- 
tional to ¢-!; it rises from zero at the plane to a maximum at the edge of 
the layer. The variation of this velocity outside the layer and in particular 
the way it falls off to zero are not determined, of course, by boundary-layer 
considerations and require a more detailed investigation of the reaction 
of the viscous layer on the flow outside. 

In the linearized solution of section 3, which is valid for small Mach 
numbers, the result in (50) is the formal expression for the pressure over 
ill the regions (i)—(iv) though it is not in a form useful for computation. 
We have seen that further study of integrals of the type (51) is required 
before this aspect can be considered complete. Nevertheless the most 
interesting physical feature is.the pressure at the wall, and that has been 
obtained in readily calculable form in (54) and is tabulated in Table 1. 
The general considerations of L.C.T. (1) regarding the form of the solution 


of equations of the type (33 


when added to our specific results give a 
reasonable appreciation of the whole field. Dissipation is important only 
in the region y = O(t), where, in the light of the pressure calculated at 
the wall, we should expect a rapid rise of pressure initially followed by 
a steady decrease to the undisturbed value. These pressure variations will 
be communicated to the flow outside this layer along the sub-characteristics 
of (33) with, in addition, a dispersive effect which will arise through viscous 
wtion. The velocity of outflow v, once p is determined everywhere, is 
determined by (20), which is effectively a diffusion equation with a known 
distribution of sources. Outside the layer in which dissipation is important, 
these changes, too, will be propagated along the sub-characteristics of (33) 
with the dispersive effect of viscosity superposed. 

So much for the actual solution of our equations. Now let us turn to the 
numerical results of section 3; their most interesting and striking feature, 
apart from the instantaneous rise, is the rapidity with which the pressure 
falls off at the plane. For : 


ir under standard conditions v, = 0-132 cm.? 
sec... a 3°3 x 104 em./sec., so that v,/a? is of the order of 10-1 sec. Referring 
to Table 1 and (56), we see that the pressure has fallen to one-half its initial 
value in about 5x 10-1 see. and to one-tenth in about 10-8 sec. 

Now the time between successive molecular collisions is of the order of 


0"! sec. Over such times the treatment of the gas as a continuum (the 
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assumption on which the classical equations we have used is based) is no 
longer possible and we must expect the first half of our pressure changes 
to require re-examination in the light of molecular considerations. Further- 
more, since our equations are based on energy considerations relaxation 
effects must also be taken into consideration (see Gunn 5). The relaxa- 
tion time for air, which is of the order of 10-* sec., represents the order of 
the time for the ‘inert’ degrees of freedom to respond to pressure changes: 
the ‘active’ degrees of freedom respond much more rapidly, in times of the 
order of 10-1 sec. We must, therefore, expect modification of our results 
up to times of the order of 10-4 sec. The early effects in the first half of 
our predicted pressure-time result are probably the most marked, though 
the others are likely to be by no means negligible. 

For times greater than 10-4 sec., our equations would be expected to be 
valid and we may notice that the excess pressures they predict are then 
negligible. Without therefore going into molecular considerations in detail 
it is reasonable to infer that the onset of the motion is marked by 
rapid pressure changes which die out in the order of 10-4 sec. Though in 
itself this inference is not conclusive evidence, it does point strongly to the 
fact that the boundary-layer solution is valid down to times of this 
order. 

The importance of molecular considerations, moreover, makes it appear 
doubtful whether any great effort should be devoted to filling in 
the gaps in the solution of the classical equations for this particular 
problem. 

It is very tempting to draw conclusions about the steady flow past a 
semi-infinite plane in the same way that Rayleigh did for the associated 
incompressible problem by replacing ¢ in our considerations by 2/U’, where 
x is distance from the leading edge. We might, in the very rapid pressure 
changes, claim to see indications of the shock wave which it seems likely 
must exist there. However, although there may be many points of similar- 
ity, this approach which, as Squire has pointed out, corresponds to lineariz- 
ing the inertia terms in the steady flow problem cannot, I feel, immediately 
admit of this extension. In the unsteady flow dissipation is the sole cause of 
outflow normal to the plate and of the pressure variations. In the steady flow, 
changes of velocity parallel to the plane also contribute to this outflow and 
in their turn also induce pressure changes. In fact an examination of the 
linearized equation for the steady flow shows evidence of some dispersion 
of pressure effects on this account though they may not (and indeed from 
experimental results it appears unlikely that they will) be sufficient to 
prevent a still very rapid change in pressure. Further examination is 
reserved to a later paper. 
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ON A PROBLEM OF INTERACTION OF PLANE WAVES 
OF FINITE AMPLITUDE INVOLVING RETARDATION 
OF SHOCK-FORMATION BY AN EXPANSION WAVE 
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SUMMARY 

The interaction, in an inviscid perfect gas with constant specific heats in th 
ratio 5/3, of a receding simple compression wave and an advancing, centred, simpli 
expansion wave is calculated by the method of characteristic coordinates, th 
interaction being taken to begin before a limit line (corresponding physically to 
the formation of a shock wave) is formed in the compression wave. The time is 
determined by which the interaction delays the formation of a limit line. It is 
found that shock formation can be delayed considerably, but not indefinitely. (For 
a formula for the time at which shock formation occurs in a simple wave, see the 


Appendix.) 


1. Introduction 
THE one-dimensional, isentropic, unsteady flow of an inviscid fluid ina 
region which is not a simple wave has been studied by several authors, 
in particular by Riemann (1), and by Love and Pidduck (2). In many 
problems of this kind shock formation could not occur; there are, however, 
others of considerable interest, in which the possibility of shock formation 
may be studied. 

The problem considered here is that of the interaction of two simple 
waves, one compressive and one expansive. The calculations are made 
for the case in which the compression wave is receding, with its front 


passing through « = 0 at time t = 0 and the expansion wave is advancing, 
being centred at x 2c at ft 0 (Fig. 1). Between the waves the gas 


is assumed to be initially uniform and at rest. 

The straight characteristics in the compression wave are converging and 
form an envelope (a limit line); the physical interpretation is that a shock 
wave must appear when or before the limit line is formed. The presence 
of the rarefaction wave delays the formation of this limit line. 

It is convenient to think of the simple waves as caused by piston 
motions. But this implies that reflection will eventually occur at the 
piston causing the compression wave and possibly also at the other piston. 
Such reflections are not relevant to the problem to be considered, and we 
do not wish to take account of them. By assuming that the velocity of the 


[Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 2 (1951)] 
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eft-hand piston is not less than the escape speed of the gas, we may ensure 
that no reflection occurs at this piston. The compression wave, on the 
ther hand, may be regarded as generated by a piston starting to move at 
time earlier than ¢ 0, from a position with x > 0, and by choosing this 
value of x sufficiently large we may postpone reflection from this piston 
for any length of time. 


\s a convenient method, however, of defining the compression wave, 





we shall specify the equation of the particle path through x = 0, t = 0. 
[his curve may be regarded as the path of a piston starting to move at 

0 at time f 0 with zero initial velocity, and if, for instance, the 
simple wave contains the particle path x }vt?, where v is a constant, 


we may refer to it as ‘of the type caused by a piston moving with uniform 





eceleration vy 
In section 2 characteristic coordinates are introduced, the dependent 
riable to be used is chosen, and the equation to be solved is found for 
perfect gas with constant specific heats in the ratio y, together with its 
lution in general terms when n 3(y+-1)/(y—1) is a positive integer. 
In section 3 the solution is expressed in terms of the boundary values, for 
ny piston motion defining the compression wave, on the assumption that 
5/3, Le. 7 2. In section 4 this solution is used to write down the 
equations of the singular lines in the characteristic plane, which corre- 


spond to limit lines in the (x,/)-plane; these equations are then specialized 


} for (i) the case when the piston defining the compression wave has uniform 
cceleration —v, and (ii) for the case when it moves under the law 
nt®, where 7 is a positive constant. In section 5 the position at 


which a limit line is first formed is considered in the case of a uniformly 
ecelerated piston. In order that it should not occur before interaction 
starts, the parameter A = 36(y—1)?Lv/az (where a, is the velocity of sound 
n the undisturbed gas) must be less than 12; the singular lines are com- 
puted (and shown in Fig. 2) for A 10. The position x, and time f¢, at 
vhich a limit line first appears (corresponding physically to the appearance 


f{ 


fa shock wave) are determined in section 5.1 for all integral values of the 





parameter A from 0 to 12 inclusive; the ratios x,/x, and t,/t, of these 
oordinates and times to the corresponding values for the compression 
wave alone, in the absence of interaction, are shown in Table 1 (section 5.1 
below) to lie between 9-7 and unity, and 3-73 and unity respectively. The 
case of a piston moving under the law x nt® is treated in section 5.2; 
the position at which the shock first appears is computed for the values 
2 and 100 of the parameter N 360L(3n)i sof. 

\ simple formula for the time at which a shock forms in a receding 
‘imple wave is established in the Appendix. 
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2. The equations of motion 
The equations governing the one-dimensional, isentropic, unsteady flow 
of inviscid gas, in terms of Eulerian coordinates x and f¢, are 
ou ou . lop 
—— ab Sf - | = \, (1) 
ot Ox pox 


( ( Ou 
= ty cP 1 pS QO, 


Cot ox Cx 


and (2) 
where u(x,t), p(x,t), and p(a,t) denote respectively the velocity, density, 
and pressure of the fluid at the point x at time ¢t. Introduce the Rieman- 


nian characteristic variables r and s, defined by 


> 
~ 


lfa l fa 

L = | —dp—}ju, (3) 

2J p 2J p 
0 

where a is the local velocity of sound and is equal to (@p/ép)', the derivative 


being taken with constant entropy. The equations become 


or. er 
ee ew le | (4) 
ot Ox 
Cs Os 
and — + (w—a) 0. (5) 
ot Ox 
Thus 
r = constant on any ‘positive’ characteristic dx = (u-+-a) dt, 
and 
s = constant on any ‘negative’ characteristic dx = (u—a) dt. 


Equations (4) and (5) may be inverted, and yield the equations 


Ox ot , 
= (u—a)—, (6) 
or cr 
CX ct m 
~ = (U+-a)—. (/) 
cs Os 
If, however, the Jacobian 
o(x.t ot ot 
J = ( —_ 2a (3) 
o(r, 8) or os 


is zero or infinite at one or more points in the region considered, the 
correspondence of the regions of the («,t)-plane and the (r,s)-plane will 
not be regular and (1, 1). 

Eliminating x we obtain the equation for ?: 


= 0: 


,, Gt _ uta) ct c(u—a) et 
2a 1 ae 


Cres or Os os or 


u and a are known functions of 7 and s and this is a linear equation for . 
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For a perfect gas with constant specific heats in the ratio y. p = Ap’, 
where A is a constant for isentropic flow, 
) al(y—1 Lu, 8 a/(y—1)—4u, (9) 


nd the equation becomes 


ot 
n (= ' ot 0. (10) 
Cres (7 S)\or Os 


where n S(y+1)/(y—1). (11) 


This equation was studied by Riemann; his well known method for the 
solution of linear partial differential equations was, in fact, developed to 
solve it. In general, the solution may be expressed in terms of integrals 
nvolving hypergeometric functions, but when n is a positive integer the 
ypergeometric functions can be expressed in finite terms and the solution 


can be written in the form 


1 (@ . @\]"[4,(r)+¢.(s) 
‘=leraletad) 


r+) 





where d, and ¢, are arbitrary functions to be determined from the boundary 
conditions. 

The complexity of the solution increases with increasing n, and it is 
easier to use the dependent variable w, introduced by Riemann, and 
defined by the equations 

cw cow 


= (uw—+-a yet. — = x (u ayt, (12) 
cl Cs 


In place of ft since Ww satisfies the equation 


C2w n L/ow . ew 
Ls 0. (13) 


Cres / Ss or cs 


This equation is of the same form as equation (10), but with n—1 in place 


+ 


ol n, so, for positive integral n, the solution is 
| n—-20 Fir)+ f(s 
Ww : : (+I) " 
(7 S$) \cr os j 


r+s 
where F and f are arbitrary functions, and these are more easily deter- 


(14) 





mined from the boundary conditions. 

In the following work an integral value of n is assumed: since we are 
mainly concerned with qualitative results, the value n = 2, y = 5/3, is 
hosen for simplicity. (A convenient method of determining f(s) and F(r) 
lor larger integral n is given by Taub in (3).) 
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3. Boundary conditions 

The interaction of the simple waves commences at t —L/a,,x = —L, 
and leads to the formation of a compound wave in the region P'0'?!' 
(Fig. 1). 





The solution in this region is of the form (14), with n 2, that is 


| , . 
w [F(r)-+f(s)]. (15) 
(7-+-8) 
and f(s) and F(r) must be determined from conditions at the wave fronts 
O'P’ and O'Q’, where x and ¢ are continuous functions of 7 and s. 
[In the compression wave the negative characteristics are straight and 


have the equation a = (u—a)t+ 7's), (16) 


where 7'(s) is a function determined by the piston motion defining the 
wave. In this wave r is a constant and equal to its value on 0’ P’, which 
1S 19 So My/(y—1). 


The equation of the positive characteristics in the expansion wave is 


x (uta)t—2L; (17) 
and in this wave s = s). Thus at the advancing front, O’ P’, we have 
[é wics r=S8o i [w : (u a Mt |, So ai T(s), 
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and at the receding front, O’Q’, 


[é w/c r|, : [a (u- a)t}, Se -2L. 
Choosing F(r,) = f(s») 0 and integrating, we have 
fis) = —(s+4) [ T(s) ds, (18) 
F(r) = 2L(s2—r?), (19) 


2L(s8—r)—(8+8) | T(s) ds 


ind thus wW So (20) 


4, Singular lines 


The correspondence between the region P’O’Q’ and the region 


S : =e 828 
0 u 59> 0 


in the (7, s)-plane, will not be regular and (1,1) if J vanishes.t When J is 
zero, either f/és or ét/ér is zero, and a line in the (2,t)-plane on which 
either vanishes is an envelope of one family of characteristics and a cusp 
ine of the other family (5). Near such a line the solution will therefore 
vield two different values of the velocity, at the same point and at the 
same time, and thus cannot represent the motion of a gas. The envelope 
of characteristics in the (x,¢)-plane is called a limit line and the corre- 
sponding curve in the (r,s)-plane is called a singular line. 


By (9) and (12), we have, for y 5/3, 
(2s r) cw/er+(s 2r) ow/os 3 ] Cow . ow 
xr t - . —< > (21) 
(7+) Z2(r+s)\or as 


(r+-s)®f"(s)—2( F’(r)+ 2f'(s))(r+s)+6( F(r)+f(s)) = 0, 


where primes denote differentiation with respect to the argument. The 
equation of the other singular line ét/ér = 0 is, similarly, 


(r+-s)?F"(r)—2(2F’(r)+f'(s))(r+s)+ 6( F(r) +-f(s)) = 0. 


The correspond e will also not be regular if J is infinite. If J is infinite either ér/éz 
must be z n this case we have a singularity in the mapping of the (2, t)-plane 

on the (7, s)-plane "he solution given above fails, though it can be modified to yield 
orrect solutior In this work, however, we assume that, for ¢ 0, the piston motion 
‘ining the compression wave is always one that has negative acceleration ; in this case zeros 


cx and @s/€x cannot occur inside P’O’Q’ (5). 
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Substituting for F(r) and its derivatives from (19), we find that the 
equations become respectively 
r( f"—4L)+r(2sf"—4f'+ 8Ls)+ s*f"—4ef'+ 6f+-12L55 = 0, (22) 
sf’ —3f—6Ls?+2Ls? 

and r= - : =; ‘ (2: 
4Ls—f 
4.1. As stated in the introduction, it is convenient to think of the 
compression wave as being caused by a piston which starts to move from 
x = 0 at time t= 0. We now specialize the formula obtained in the 
previous paragraph for two particular piston motions. Suppose firstly 
that the equation of the piston curve is x = —}vf*, where v is a positive 
constant. Let «,,(s) and t,,(s) be the values of 2 and ¢ at the point at which 

a typical straight characteristic, s = constant, meets the piston curve. 


Then 
—vt, u S—S8p t, = a Ly = — yy 
Substitute in (16) to obtain 
T(s) = (5s—6ss,+-s?)/6v; (24) 


and with this value of 7'(s), from (18), 
f(s) 
Substitute for f(s) in (22) and let 


A = 36Ly/s?. (26) 





(5s4 — 4333, — 6s? + 4ss3 - s) 18p. (25) 


The equation (22) of the singular line, é¢/és = 0, corresponding to an 
envelope of the negative characteristics, is then found to be 
r2(60s?— 2488, ee 128° iL. 2As?) }. r(40s* 4. 2488, — 1688 _ 4X83?) 1 


+ 10s4+ 885+ 6s§—6As§ = 0. (27) 


By substituting for f(s) in (23), we find that the equation of the other 
singular line is 

- 594 6385 —Assp 8ss3 + 3As8 — 384 (28) 

2083+ 2Ass2— 12s?s,— 1282+ 483 * . 


Secondly, we obtain the solution in the case when the piston law is 


a 7t3, where 7 is a positive constant. By the same method as used 

above, we find that ion (s—s,)¥(30- S0). si 
(27)! 

and f(s) = —4K(9s?+ 1088) + 82)(s—sp)!, (30) 

where a (31) 


~ 90(3n)* 
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5. The position of shock formation 
The physical interpretation of the appearance of a limit line is that the 
sentropic flow breaks down and a shock is formed. The best estimate, in 
the present stage of knowledge, for the position and time at which the 
shock wave is formed is the position and time at which a limit line first 
ppears. We therefore calculate the lowest value of ¢ on the limit line. 
The singular lines given by (27) and (28), for A = 10 are shown in Fig. 2. 





Fic. 2. 


The singular line ABC is the line ét/és = 0, and CFD is the line 
tor = 0. In region I, ét/@s and ét/ér have the same sign as in the simple 
vaves: that is. et/er 0, et/es 0. Moreover, é7t/erés ~ 0 (except at B); 


ence in 


Region I] ot/d6s <0, ot/or < 0, 
Region III: cet/eés < 0, et/er > 0, 
Region IV: cét/és > 0, et/ér > 0. 


herefore, the time has its minimum value on CFD, and its maximum 
value on ABC, at B. On the singular line BGC, which corresponds to an 


envelope of negative characteristics in the (2,t)-plane, we have ét/és = 0, 


‘cr > 0, and we find that ¢-» —oo as we approach C along this line. 
Hence the earliest appearance of a limit line is at ¢ 2, but clearly 


his can have no physical significance, since interaction did not begin until 

Lia). We must restrict ourselves to that part of the (x,t)-plane in 
vhich our isentropic solution can have physical significance, that is the 
region which can be reached from the simple waves without crossing a 
singularity. This region corresponds to the region in the (r,s)-plane 


ounded by the line 7 = s», a segment of the line s = sy, and the singular 


ines AB and BD. Shock formation will occur at the lowest value of ¢ on 


092.14 


N 
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the limit lines corresponding to A B and BD, i.e. at the point corresponding 
to A. The isentropic solution which we have found in section 4.1 will hold 
only in the region of the (r, s)-plane bounded by OA, r = 85, and the curve 
through A which corresponds to the shock path. 


5.1. We now investigate the position of shock formation for com- 
pression waves of the type caused by a uniformly accelerated piston, for 
general A. It is shown in the Appendix that, in the absence of interaction, 
a shock would form at the head of the compression wave at x 12). 
in order that the shock should not form until after the interaction has 
commenced we must choose 

L<4s/v, that is, A < 12. 

To determine the position of shock formation we have to find the lowest 
value of ¢ on the singular line ét/eés = 0, as in the case A 10. 

If we take s as a parameter along the singular line, and let the equation 
of the singular line be r R(s), we have, differentiating ¢ along the singular 


ume, dt : (él) dR [ét\ dR 
ds Fal = ( | ds ( ‘) ds~ 

The zero of et/cr can be shown to correspond to a maximum of ¢ on the 
singular line, as in the case A 10, and to occur at an end-point of the 
singular line segment under consideration; let this point be B(A). A mini- 
mum of ¢ on the line can occur only at a point where dR/ds = 0; if there 
is no such point the lowest value of ¢ will occur at s = Sp. 

We have found numerically, for A = 10, that dR/ds < 0 on the segment 
AB of the singular line. Suppose a zero of dR/ds occurs for A = A* inside 
the corresponding segment, but for no value of A between A* and 10. Then, 
dR/ds being a continuous function of A (except at B(A)), dR/ds < 0 on the 
segment; moreover, dR/ds is a differentiable function of s (except at B(A)), 
and hence d?R/ds? also vanishes where dR/ds = 0. 

From (27) we find that dR/ds d? R/ds? 0 implies R s: but 
r s is the equation of the vacuum line, and hence there is no point 
on the segment of the singular line at which d?R/ds? = dR/ds = 0. Thus 
as we increase or decrease from the value A 10 a zero of dR ds must 
first occur at s = Sp. 

Differentiation of (27) with respect to s yields 
R?(120s—24s,)+ R(120s?-+ 2482 — 4s?) + 4083 + 883 + 

dR 


7, (2:R(60s?— 24889 1252 + 2As?) + 4053+ 24552 — 4382 — 1653)} = 0 


whence, 


dR 24 R?(8,) +89 R(s9)(36—A)+ 1288 
(Ts). " 8o{ R(sp)(12+A)+12—At* 
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dA— 12 
Also, from (22) R(s,) aT 8) or Sp: (33) 
Substituting the first of these two values in (32), we find that 
(dR/ds), ma 0 


A\? 26/A' 1—0 
aun (73) 3 {73} : ; 
There is one root of this equation in the range (0,12), namely 1-404, and 
for this value of A, dR/ds = 0 on the limit line at s = 8, i.e. A* 1-404. 

Thus, if A lies in the range 12 >A 1-404 shock formation occurs at 
the head of the compound wave; if 1-404 > A > 0 shock formation occurs 
nside the wave. 

The values, r, and s,, for which ¢ has its lowest value on the singular 
ine segment considered are given in Table 1 for integral values of A. For 
1-404 < dA ae s, and r, is given by (33); for A 0 and A 1, the 
singular line has been computed to determine the zero of dR/ds. The 

responding values, x, and ¢,, of x and t, which provide an estimate for 
the position and time at which a shock is formed, are also given in Table 1; 

and ¢, denote the position and time at which a shock would be formed 
f the rarefaction wave were absent (cf. Appendix). The values of x, and 
ire found from (20), (21), and (25) by substituting r r, @= 6 ie 

$()4 A 12 


L(3p?—4y+6)/8, t, = 3L(1—p?)/8s,, 
where pu ZIA 
TABLE | 
3°73 1 7 
7 122 7°405 
$2 f124 
1°562 2-688 
1*204 I ) 
1125 1°375 
1-074 1°24 
1-042 1°125 
I 21 I ( 
I I 5 
I I 18) 
I I > 
rhe case A = 0 is a limiting case; A may tend to zero in two ways. If 


L-+( and » remains finite x. and ‘. tend to finite values. thus we cannot 
ostpone shock formation indefinitely if the piston has non-zero accelera- 
tion; if y+ 0 and Z remains finite x.—> © andt > oO. 
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The values in Table 1 suggest that, for the values of A for which the 
shock forms at the head of the wave, the ratio t,(7,—~2,)/x,(t,—t,) is inde- 
pendent of A, and is equal to 3. It can be shown that this result is true 
for any compression wave if (a) y = 5/3, (b) the shock is formed at the 
head of the compound wave, and (c) if in the absence of interaction the 
shock would form at the head of the compression wave. 

5.2. The position of shock formation when the compression wave is such 
that, in the absence of interaction, shock formation would occur inside the 
wave has also been studied in one special case. The case considered is that 
corresponding to a piston moving under the law x nt?, 1 > 0. The 


parameter corresponding to A is N 360L(3y)'sg ?, and to ensure that 
shock formation does not occur before interaction it is sufficient to take 
N < 217. In this case, in the absence of interaction shock formation will 
occur on 8 11s,/9. 

The equation of the singular line in the (7, s)-plane is found by substitut- 
ing the expression for f(s) from equation (30) in (22). This line has been 
computed for the values V = 200 and N = 100 and the zero of dR ds 


determined. The results are given in Table 2. 


TABLE 2 


\ / 
\ SQ V./ Sq l/l, XQ/X¢ 
100 I°2O01 | O°315 I*10o 1*274 
} 


200 I°215 | 0-850 1018 | 1043 

The author wishes to thank Professor 8S. Goldstein for suggesting this 
problem, and both Professor Goldstein and Dr. R. E. Meyer for their 
helpful criticism of this paper. 


APPENDIX 
Shock Formation in a Simple Wave 


It is shown in (4) that when a piston moves into a gas at rest there is a limit point 
on the front of the wave generated by the piston if and only if the initial acceleration 
of the piston is non-zero. This need not, however, be the first limit point to appeal 
during the motion, and hence need not provide an estimate for the time and position 
at which a shock wave first appears. A simple formula for the time at which a limit 
point first appears may be established as follows. 

Let the speed of sound in the gas at rest be a, and let a receding wave move into It. 
Then r = 8 = ao/(y—1) by (9), and the equation of the straight negative charac- 
ee x (w—a)t+ T(s). (34 
The limit line is an envelope of the characteristics, and hence it is given by (34) to- 
gether with dx/ds dt/ds 0, i.e. by (34) and 


t 2T’(s)/(y+-1). (35 
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If the motion of the piston (or indeed, any particle path) is given by 2 = z,(s) 


- 
id ft t,(s), then mm 
and T'(s) x,(8)—(u—a)t,(s), 
ind dx,,/dt u, the speed of the piston. Moreover, by (9) since r Bo, 
dt,,/ds dt ,/du 1/b,, 

where 6,, is the acceleration of the piston. Hence 

T’'(s) (y at 1)t,, 2 D,, 
ind, by (35), the time at which any one of the straight characteristics s constant 


meets the limit line is 


t t,—2a/{(y 1)b,}. 
The earliest time, t,, at which a limit point appears is therefore equal to the lowest 


ilue of 


the piston curve (where ¢ is the time, a the speed of sound, and b the acceleration 
f the piston), and this limit point lies on the straight characteristic through the 


piston curve where the lowest value occurs. If t—2a/{(y+ 1)b} is a decreasing function 


fton the piston curve at the point where the piston starts its motion, the shock 
| not form at the front of the wave, even if the initial acceleration is non-zero. 
Similarly, the time at which a limit point first appears in an advancing, simple 


pression wave is equal to che lowest value of t+ 2a/{(y-+ 1)b} on the piston curve. 


When 2 svt? and y 5/3, 7's) is given by (24), and by (35), 
t 489/v. 


From (34), the corresponding value of x is found to be 


t $85 /v. 
When nt and vy 5/3, T(s) is given by (29), and by (35), 
t (37) 4(9s 78, )(8 89) } 
s has a minimum when s lls,/9, and 


7 } 
(49/9)? 
} 


23(a3/y)?/36. 
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ON THE IMPULSIVE MOTION OF A FLAT 
PLATE IN A VISCOUS FLUID 


By K. STEWARTSON 
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[Received 21 February 1950; amended 23 November 1950] 


SUMMARY 
The motion in the boundary layer which arises when a semi-infinite flat plate is 
impulsively started from rest with velocity U is investigated. It is found that the 


velocity field is independent of x, the distance along the plate if Ut < x, and at 
Ut = x the flow has an essential singularity, thereafter depending on x as well as t. 
As t-> © the influence of ¢ is shown to die out exponentially. Solutions are also 


obtained using two approximate methods and they both agree in requiring the 
velocity to be independent of x for a finite time. 

The same problem is also solved for a compressible fluid. It is found that in this 
case the outflow is increased and this may lead to the formation of two compression 


waves. 


1. Introduction 

THIS paper is concerned with the motion of the fluid in the boundary layer 
associated with a semi-infinite plate, which starts to move impulsively 
along its length with constant velocity U’ in a fluid at rest at infinity. 
Previous papers on this subject, of which Goldstein has given a connected 
account (1), consider the impulsive motion of cylinders for which the 
relative velocity U, of the fluid outside the boundary layer is not constant, 
and in which separation occurs fairly quickly, before the boundary layer 
has had time to settle down to its steady state. The solution is obtained 
in the following way. Initially the boundary layer has zero thickness and 
therefore at the beginning of the motion the diffusion is far greater than 
the convection and the influence of the pressure gradient. Consequentl) 
the velocity in the boundary layer is a function of y, the normal distance 
from the surface of the cylinder, and the time ¢, multiplied by the relative 
main stream velocity, calculated from the inviscid theory. A second 
approximation is then obtained by substituting it back into the boundary- 
layer equations. The process may be repeated as many times as desired, 
the only restriction being that the boundary-layer equations are invalid 
beyond a point of separation. A peculiar feature of the method is that if 
as in the problem to be discussed in this paper, U is constant, all the terms 
after the first are zero, so that the leading edge of the plate apparentl) 
exerts no influence on the motion in the boundary layer. Such a conclusion 
is paradoxical since the flow eventually steadies down to the Blasius 


(Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 2 (1951)] 





solut 
in th 

Be 
thro 





simp 
linea 
wher 
veloc 
leadi 
It su 
of th 
if Ut 

TI 
a sol 
if Ut 
but | 
skin 
metl 
chan 
the | 
com 

Tl 
unfo 
toa 
mat 
dent 
esse 
and 
abor 
pass 
ism 
is $I] 
solu 
this 
that 
may 
argy 
find 
only 
esse 


an ¢ 








SIVé 


finity 





econ 
dar 


sired 


nvail 
hat 
terms 
rent } 


‘Jusiol 











ON THE IMPULSIVE MOTION OF A FLAT PLATE 183 


solution. It is of interest, therefore, to examine more closely the solution 
n this case. 

Before examining the full boundary-layer equations, some light may be 
thrown on the general nature of the flow in the boundary layer by two 
simple approximate methods. The first, known as Rayleigh’s method, 
inearizes the boundary-layer equations and is only valid at the outer edge 
where the fluid viscosity is approximately the same as the main stream 
velocity. The solution it gives is independent of x, the distance from the 
eading edge along the plate, if Ut < x, and is independent of t if Ut > x. 
It suggests that the fluid velocity at (x, y) is unaffected by the leading edge 


of the plate if Ut < a, and is independent of the time of the initial pulse 
fUt>z 

The other method, which makes use of the momentum integral, gives 
, solution which is independent of 2 if Ut < 2-65x, and is independent of t 
{Ut > 2-65x. Both these methods inevitably over-simplify the problem, 


but in view of the success of the momentum integral in determining the 
skin friction in accelerated flows, and of the validity of the Rayleigh 
method at the outer edge of the boundary layer, it seems likely that the 
change over from ¢ to x as a principal independent variable will begin at 
the outer edge of the boundary layer when Ut = x and will be nearly 
omplete when Ut = 2-652, having spread down to the plate. 

The boundary-layer equations for this problem are complicated and, 
unfortunately, the conclusions derived from an examination of them are 
toacertain extent tentative. Since the straightforward method of approxi- 
mation breaks down, it might be expected that the flow would be indepen- 
dent of x for a finite time and that then 2 would enter by way of an 
essential singularity. We find that a solution of this form can be found, 
ind that in it the critical value of ¢ is x/u. The velocity in this solution 
ibove any point of the plate is therefore independent of x until that point 
passes through the original position of the leading edge. The influence of x 
is most marked near the outer edge of the boundary layer when (Ut—2x)/x 
is small and increases inwards with (Ut—a)/z. When Ut/x is large the 
solution as a power series fails in the same way as when Ut/x is small. In 
this case a solution is found with an essential singularity at 2/Ut 0, so 
that the flow would eventually steady down to the Blasius solution. It 
may well be that there are many other possible solutions, for none of the 
rguments used vbsolutely preclude them, but it has not been possible to 
find any excep those with essential singularities at Ut = x and 2/Ut 0 
only. There are actually an infinite number of possible solutions with 
‘sential singularities at x/ Ut = 0, but with one exception they all contain 


n oscillatory component. 
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The solution we find does agree with Rayleigh’s method if Ut < x, as 
indeed we should expect, for the latter is then an exact solution. It is only 
when Ut > x that the flow becomes mathematically complicated, but stil 
the general features are as outlined in the last paragraph. 

The discussion above has tacitly assumed that the fluid was incompres- 
sible, but a solution of the same form may also be obtained when com- 
pressibility effects are taken into account. If we transform y by means of 


y 
y_ [ey 
A Po 
0 


(1.1) 


where p, py are the densities at x, y, t, and infinity respectively, and assume 
that the viscosity is proportional to the absolute temperature, the stream 
function can be shown to satisfy the same equation whether the fluid is 
compressible or not. The solutions are slightly different, however, since 
the compressible solution is completed by transforming back from Y to y. 
In addition tothe outflow when Ut > x, which occurs in the incompressible 
solution, there is an outflow set up at t = 0. 


2. The statement of the problem 

The problem may be regarded as one in which a velocity is impulsively 
set up in the main stream, the plate remaining at rest, and we shall treat 
it in this way here. Let us take as origin a point on the leading edge of the 
plate, the z-axis along its length, the y-axis perpendicular to it, and let the 
motion begin at ¢ = 0. If uw and v are the components of fluid velocity 
along and perpendicular to the plate, then for an incompressible fluid the 
boundary-layer equations are 


ou . ov 
|. — 0, (2.1) 
cx OY 
Gu eu 1 Gp nu 
and 6 — 49 ae —— Sy (2.2) 
Ox oy p Cx oy 


When the fluid is compressible the equations are slightly different, and 
they are given in the last section (equations 7.3—7.6). 
The velocity in the main stream is constant and equal to U’ relative to 
the plate, so that @p,/éx = 0 and the boundary conditions are 
u=v=0, wheny=0,r>0,t>0; 
woe 7, when y > 0 and either t > 0, x = 0, ora >0,t = 9%; 
u—> U, asy>o,r>0,t>0. 
The two approximate solutions of this problem will now be given before 
considering the boundary-layer equations in detail. 
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ON THE IMPULSIVE MOTION OF A FLAT PLATE 


3. Rayleigh’s method 

The basic assumption of this method is that the variation of velocity 
from the main stream value is small, so that the equations may be linear- 
ized. Such an approximation is valid only at the extreme edge of the 
boundary layer, since as we move into the plate the velocity tends to zero. 
Nevertheless the method has provided some useful qualitative information 
about the nature of the flow in the boundary layer. The boundary-layer 


equation then becomes ou oly 
cl 7 C v ‘ 
LU - =v—, (3.1) 
ot ox oy? 
neglecting products of v and (U—u). 
Write f = y/,/(vt), 7 Ut/x, and then from a dimensional argument u 


is a function of ¢, t only. It will therefore satisfy 


o*u C ou o, OU ‘ 
=+=— (r—r?) —, (3.2) 
OL" 2 OT 
with boundary conditions u = 0, when { = 0, and u > U, as > © or as 


7 ©. It may be shown that there are no other analytic solutions of the 
equation. Let us search for a solution of the form 


x 


u = | f(s,z)sin fs ds. (3.3) 


0 
There will be a solution of a similar form with sin sf replaced by cos s¢, 
but this will vanish identically, since « = 0 when ¢ = 0. Substituting in 
the differential equation, we find in the usual way that equation (3.3) satis- 


fies equation (3.2) if 


sf sin s¢], 0 and 3s J + 4f+(r—r?) gy | s?*f = 0. (3.4) 
és Or 
Hence u (= (- sin sl ds, 


where g is an arbitrary function. As t > 00, u—> U and hence 
g(+s*) . (3.5) 


As (-+ «©, u-> U, which is satisfied unless 7 > 0 simultaneously, so that 


7@ is constant. It is still satisfied even then if g(—s?) = 2U/z. Therefore 
if 7 l, 


OU f4 5 Gs sale —* 
u sin sf ds U erf : = Uerf y ‘ (3.6) 
7 8 2 2, (vt) 
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and ifr > 1 


2U fer. , tt 5 U 
u : sins€ ds = U erf mf eff / 3.7) 


7 8s Z \2 va} , 


The significance of this result is more readily seen if we refer the solution 
to axes fixed in space, and coinciding instantaneously at t = 0 with those 
we have already chosen and moving with the plate. 

Distinguishing the new axes by the suffix 1, the solution is 


where A = tifx, > 0 and A = t+2,/U if0 >a, > —Ut. 

The fluid above that part of the plate. to the right of the initial position 
of the leading edge is behaving as though the fluid were of infinite extent 
and impulsively started from rest at ¢ 0, while that to the left is behaving 
as though the plate were finite and had always been moving with uniform 
velocity. This solution satisfies the full equations of motion so long as w is 
independent of x, and hence may be expected to be correct then. Confirma- 
tion of this view will be found in section 5 below. One drawback to this 
method of approximation is that it cannot be improved directly. If the 
flow were steady the next approximation would be less accurate than this 
one, and in our problem the next approximation is discontinuous at Ut = x. 
Such a result is not surprising if the essential singularity we find in section 5 
does actually occur. The drag at any point of the plate may be deduced 
from the value of éu/éy when y = 0. According to this approximation 


Cu U , x 
: 2 where B min{f, —)}. (3.9) 
cy \(7vB) l 
Equation (3.9) is exact if Ut < x, but for large t it is known to over- 
| £ 


estimate the drag. Further light on the drag variation may be thrown 
by the momentum integral. 


4. The momentum integral 


Using Goldstein’s notation (2) we define 


; 
= | r(! 7) dy; oF | (1 = r) dy, (4.1) 
0 


where 6 is the thickness of the boundary layer. 


0 


The momentum integral is then 


p72 oF 1U 08, 
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If now we assume with Lamb (6) that 
ne ee ; 
3.7 u U sin— (0O<y< 9d), (4.3) 
20 ‘ 
a ; ae eu) aU 
tion then 0; 0-56050 Uv 0-1370, =? (4.4) 
oy y-0 <0 
hose . 1: 
and the momentum integral 1s 
Farol = FT 082 - 
0-363 0-137U - TV. (4.5) 
cl Ox 
(9.8) | This is a great simplification of the problem, since we assume that the 
| velocity profile is determined by a scale factor only. The boundary con- 
; ; ; ; 
ditions are that 6 0 when ¢ 0 for « > 0, and when x 0 fort > 0. 
Ol The solution of this equation is 
ter +e 77 mY 
: 5 2-94, (vt) if Ut 2-652 | 
Ving > 
var eel ; (4.6) 
forn 5 1-79 (75) if Vit > 2-65ar | 
wis N 
‘ma The corresponding formulae for the drag are that 
this cu . v 
+ (p 0-534pl if Ut wane 
ULI \ cy AN t ai 
this . i ‘ (4.7) 
* cou ’ V °° TT - 
: [p 0-328pl -. 2 > 2-65 | 
; CY 0 N hor 
ond rm . 
es These results, unlike those obtained by Rayleigh’s method, give good 
. greement both with Rayleigh’s solution for the infinite plate (¢ + 0), and 
; 
with Blasius’s for the flat plate (f + 00). The correct numerical values are 
3.9 ):565 and 0-332 respectively. The method is also known to give good 
agreement with steady problems in which separation does not occur, and 
vel this taken in conjunction with the agreement for small and large t suggests 
rown | that equation (4.7) is a good approximation to the truth, although the 
transition is certainly not so abrupt. 

Both of these approximate methods require t to go over to x as indepen- 
dent variable at some value of 2/Ut, which is an over-simplification of the 
transition from impulsive motion to steady motion. In the next two 
sections we shall use the full equations of motion to investigate more 

(4.1 closely the way in which this transition takes place. 
9. A solution in the neighbourhood of Ut = x 

When ¢ is small the effect of diffusion far outweighs the effect of con- 

vection and the influence of the pressure gradient, and therefore, in par- 
ticular, the effect of the leading edge. Hence the fluid will move according 
(4.2 4 > ‘ ; > : . : —_— 

0 Kayleigh’s solution for the impulsive motion of an infinite flat plate (4). 
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Since we cannot approximate by means of a power series to the solution 

for larger t, there must be an essential singularity in the solution for some 

value of ¢, and in this section we propose to investigate its existence. 
Let us take as independent variables 


‘s y Ut 
Peay eel Pa os (5.1) 
(vt) x 
and write u U,/(vt) dl v —U,/(vt) “dl (5.2) 
ey ox 


On dimensional grounds uw must be a function of { and 7 only, and there- 
fore, substituting in the equations of motion (2.1), (2.2), & satisfies 


Sy Lap ( ped 4. Bab Orb 
= = |T—T* 4+ 7% 


Uy 


(5.3) 


bo 


of oc* 06] 0COr © 00? er’ 
with boundary conditions 
2 = ob/eC = O when €=0; a@/el>lasf>o, orast>o. (5.4) 


As t > 0, 7 + 0, and ys > &, independent of 7, and is given by 


ous m4 1 
—~ =: orf—. (5.5) 
rar 2 
: nae Cus 2 exp(— (2/4 tins 
When € is large =e ie § P\ ; ) (5.6) 
: (ae CQvVT 


8 


and has an essential singularity at ¢ 

Equation (5.5) satisfies the differential equation (5.3) and the boundary 
conditions, except that it does not tend to unity as 7 tends to infinity. 
The obvious next step is to assume an expansion for y% as a power series 
in t whose coefficients are functions of ¢, and whose leading term is yp. 
However, on substituting in the differential equation it is found that all 
terms vanish except the leading one, which is independent of 7. This 
implies that the solution must have an essential singularity at some value 
of r. Such a singularity is not necessarily impossible physically and indeed 
may be highly desirable, if all derivatives with respect to 7 vanish there. 
One example is given by ¢ = 0, 7 < 0, and % = e-"7, 7 > 0, which has a 
smooth transition at = 0. Little is known about the occurrence of 
essential singularities in differential equations, but they are only expected 
when the coefficient of a leading term of the differential equation vanishes. 
Consider the differential equation in the neighbourhood of such a singu- 


larity. If we write ob = the(L) +H, (7, 2) (5.7) 


and neglect 47 we have 
Oy, ¢ C7, (> . ae 
7 ‘ 
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The leading terms of this differential equation are @y,/6¢? and é7y5,/é7eC, 
and the coefficient of the second vanishes when + = 7? éy,/0, i.e. r = 0 
ort > 1. However, when + = 0 the coefficients of the terms arising from 
the finiteness of the plate are of order 7? (these terms come from the 
differentiation with respect to x in equation (2.2)), while the others are of 
order rt at most. Hence these terms exert no influence on the first-order 
solution and it remains independent of r. 

Let us now consider 0 7 < 1, so that the coefficients do not vanish. 
In this case there is no general proof that an essential singularity cannot 


occur, but a simpler equation may be used to show that it is unlikely. 


Conside , ea ; 
on ider o*dh Cod ] Ch (5.9) 
_ : o.% 
Oe 20 Zor 
where ¢ vanishes when ¢ = 0, when + x, and when ¢ 2. This equation 


may be obtained from (5.8) by neglecting ¢y,/é7 and simplifying the coeffi- 
cient of é*s, eer. The general solution of this equation for 7 > a may be 
obtained as in Rayleigh’s method (section 3) and is 

eT—%d, | sin Ct f(te*-7)e dt | cos Ctg(te*-7)e dt, (5.10) 

0 0 

where f and g are only restricted in so far as the integrals are finite. Then 
since d = 0 when € = 0, g = 0, and since ¢ = 0 when 7 = a, f = 0 and 
therefore 4 0. 


Having shown that it is unlikely that an essential singularity will occur 


before s = 1, we shall now investigate the possibility of such a solution 
at 7 = 1. It must be emphasized at the outset that the results are more 


tentative than conclusive. On the other hand, it seems likely that the only 
possible finite value of 7 at which an essential singularity can occur is 7 = 1. 


The difficulty about this equation near 7 1 is that when € is large and 
T 1, the coefficient of 
O2/, ) (2) 
Y1 nid ; : 
~ ——exp| —=], 
CQCT ONT 4 
irom equation (5.6) and hence has also an essential singularity at € = ©. 


We shall now investigate the behaviour of such a solution as 


(7 1) (1 3) > CO, 
Cg 


and then link it up with a solution for more moderate values of 


, n/t *) 
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7—1 being invariably small. Let ¢ — 00 in equation (5.8), and keep 7—] 
small. The dominating terms of the equation satisfy 


ay 
a 1) - Pr (5.11) 
oc OTe’ 


9 


C su, C Cp, 
1° 


um 


C 
Hence for ¢ large the solutions of (5.8) of order (7—1) are dominated by 
(r—1)e-$"/4, (r—1)/C, and 7r—1, (5.12) 


while those of order (r—1)? are dominated by 


(r—1)2f2e-$"/4, (r—1)?/Z3, and (r—1)?. (5.13) 


Consequently any solutions we obtain of order (r—1) or (r—1)* must be 
linear combinations of these as ¢ > oo, and hence their derivatives with 
respect to ¢ must tend to zero as €-> 00 as required for %,. Since the 
derivatives with respect to 7 must also vanish at +t = 1 (so that v will be 
continuous), we must choose our solution carefully so that the third term 
of equation (5.12) cannot occur. The behaviour of %, for moderate values 
of (r—1)/(1 
ducing a new variable instead of ¢, where 


6ys,/€C) with s—1 still small may now be examined by intro- 


8 > — 
b —-e-* - (5.14) 
3A a 


um 


Substituting in the differential equation and retaining only the dominating 
terms for 6 small, we find that uy, satisfies 
3 2,/ } : i 07y - ] 1 
a OW 9 OW ) ou Oo-w ous ss Ses 
53 1 2b2 i L b | b 1 (5.15) 


cb cb? 2logb eb éreb\ — logb 


CT 
The boundary conditions are that 4, — 0 when 7 > 1+ and also as b > ~; 
xs, /0f > 0 as b> @ and yf, is asymptotic to a linear combination of the 
solutions in equations (5.12) and (5.13) as 00, ie. b> 0. A major 
difficulty is the presence of log b in two of the coefficients of the differential 
equation (5.15). Since log) varies more slowly than b, some information 
may be extracted by setting logb = c, a constant, in these terms. It is 
precisely this assumption which makes our results rather tentative. Any 
dominating terms that are found may have to be multiplied by some power 
of logb, but this should not affect the general conclusions. With this 


assumption a solution may be found in the form 
ob, (r 1)?f(s), (5.16) 


where s b/(r—1) and f satisfies 


asf nar 2) of 9 
<< C  , = c ; 2s , in at 
3 J { J (20% |_ 3 _| J 252 a 2sf 0. (0. 17) 
es 3 = 
ds* és? c Os c 


The form of equation (5.16) has been used for y¢, to eliminate the possibility 
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ON THE IMPULSIVE MOTION OF A FLAT PLATE 19] 
that ys, ~ t—1 as > 00 which implies that v would be discontinuous at 
r 1. Near s 0, f is a linear combination of 

s(—logs)?, (—logs)!, and 1 (5.18) 


from equations (5.12) and (5.13). All of these satisfy the boundary con- 
ditions on ys, as b and therefore s > 0. On the other hand, when s is large 


the dominating part of equation (5.17) is 


53 } 53 £ ui 22 I 1 sf = 0, (5.19) 
Cs? Cs" Os ; 
the solutions of which are / 8, f s? and 
e-§ (r—1)4 b - 
j~ > or yn = exp : (5.20) 
s* b? (r—1) 


Of these only the latter is acceptable as s — , since f tends to zero expo- 
nentially there. Hence we have found a solution satisfying all our con- 
ditions; for as tr > 1+- for fixed f, s > 00 and f—> 0; as +> 0, s > o and 
again f > 0, so that wu, v are zero at the wall; as ¢ + oo for fixed 7, s > 0 
and all of the three possibilities (5.18) satisfy the conditions that éb,/ef > 0 
ind js, er be continuous there. Its range of validity may now be improved 
by repeated substitution into the differential equation and the arbitrary 
constant may be determined by the boundary condition as 7 > 0. 

Although we have found a solution with an essential singularity at 7 = 1, 
there is still a possibility that an essential singularity may not be developed 
until Ut vv, Where a |, and that % is independent of x if Ut/x < a. 
We may only say that in this case the methods which have proved success- 
ful in the rest of this paper lead to a contradiction. Hence although the 
possibility cannot be ruled out, its occurrence seems unlikely. 


6. The nature of the solution at large times 


Let us take as independent variables 


U x 
é y | ) n ~ (6.1) 
A) \vx Ut 


C " . 
nd write u Pr (1 rl o v - py (val yt. (6.2) 
cy CXL 
On dimensional grounds wu must be a function of € and Conly, and therefore, 
substituting in the equations of motion (2.1), (2.2), d satisfies 
7 1)~ (6.3) 


3 » ¢2 i “ 2 
2 2 o&- o& ¢ 


( sd, d c 24 c db 9 oh 9 7d od 
T T ae / 2; 
c&On o&" on 
with boundary conditions ¢ = é¢/é€ = 0 when € = 0; and é¢/é€ + 1 both 


S¢— cand as 7» >. As n > 0, t + © and the motion becomes steady, 
in Which case the correct solution is Blasius’s, for which ¢ is a function of € 
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only. Denoting this solution by ¢o, then ¢po satisfies all the boundary con- 
ditions except that it does not tend to unity as 7 tends to infinity. Again, 
it is not possible to expand ¢ as a power series in » whose coefficients are 
functions of €, since then the term arising from the time variation of x js 
of smaller order than that of the others. Again this suggests that we look for 
an essential singularity, and the work of the previous section suggests that 
we look for one at 7 = 0, since it is unlikely to arise at any other value. 
[t must be borne in mind in this section that we are now trying to introduce 
a time variation into the solution, whereas previously we had been trying 
to introduce a variation in x. Such a singularity is likely at » = 0 for the 
coefficient of €76/@&6n, n(é¢/6€) — y” has a simple zero at » = 0, except that 
if € = 0 as well, it has a double zero. 

Let us consider a possible form for such a singularity at 7 = 0, bearing 
in mind that, as in all boundary-layer problems, we cannot be absolutely 
certain that our solution will occur in practice, since no uniqueness theorem 
has been proved. Indeed there are an infinite number of possibilities even 
for the form we find. However, we can be much more confident of the form 
of such asolution than we could be in the previous case, because we shall not 
have had to assume that any slowly varying function, such as log }, is con- 
stant, and we have been able to obtain an explicit form for the dominating 
term. 

Let d = do(€) +f(E, n); (6.4) 
where f+>0 as n> 0 and as +0; af/e& > 0 as € +0 and as >. 
Further introduce a new variable @ instead of € defined by € = 6y. Then 
if 7 is small f satisfies 
do 
n*00? 


Ody \( &f lof O@f dbo of @ of 7 
a — ee aq)» (8.9) 
o€ c0én = 00 7» OP? 0&7 \ ° an o6 


St Fe Oe Me 


7? aps | 270 12002 | 


neglecting squares of f. Consider the following function 


, ] od a’ = 
(€,) = — —°9(0, n)exp| ——-], (6.6) 
ME 9 @ o& Sait 978 
which has an essential singularity at » = 0, ifa > 0. It will appear later 
that as @ > 00, g ~ 6. Hence although when 7 is small and @ is moderate 


Lebo 
0 a€ 


we still need this term, because without it f would not tend to a finite limit 
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as £ > 0. For moderate values of 0, since 7 is small, é¢,/e€ = 0-332€, and 


therefore g satisfies 


3, a> . Aa? F 
o'”d 1 RY 1) 9 : ; g- O(7?), (6.7) 
oH? 5, | C p § 


where A = 0:332. From the boundary conditions on f, ég/¢@? must also 


vanish when @ 0. The appropriate solution is 
2/ » 3 
o-”9 see “}'| (6.8) 
A 3 


cb? 


x (AO—1)!K, 


v0 





where K, denotes a Bessel function of order 4. This form of the solution 
is also valid when A@ |, because it is merely a power series, with integral 
indices, in A9—1, but it may also be expressed in terms of the Bessel func- 
tions J, and J_,, (3). We may now determine a from the condition that 
é9/e0® — 0 when @ = 0. For expressing é*g/06* in terms of J, and J_,, if 
9/c8 = 0 at 6 = 0, then 


2ai 

lytJ.(y)+ytJ_.(y) 0, where y = ——-. (6.9) 
dy F | 3AV3 

Since these two functions are oscillatory when y > 0, there are an infinite 
number of roots of this equation of which the least is y = 0-613 or 
a = 0-648. If we considered any other root, g would be oscillatory when 


\9 < 1, which is not very likely to occur. When @ is large 


3 
a(S): (6.10) 


g and ég/e@ may be obtained by integration and using the boundary con- 


od l 
~ exp 


C A2 H 





ditions when @ = 0. From equations (6.10) and (6.7) it follows that 


~AA +o(1), (6.11) 


ct 


where A is a function of » only, and 


g ~ (AQ—1)A+0(I). (6.12) 
Hence for moderate values of €, and » very small, 
od as : 
f~A exp seca | (6.13) 
c& 97% j 


and this may be verified to satisfy the differential equation, and tend to 
a finite limit as € > oo so that the boundary condition there is satisfied. 
We have therefore obtained a possible form for the solution near 7 = 0 
with an essential singularity at » — 0. The boundary condition as y > « 
may presumably be used to determine A. According to this solution the 


influence of » on the velocity parallel to the plate is greatest near the plate, 


dying out like é%4,/eé? as € > o0 and like exp(—0-03/n3) as n > 0. 


5092.14 


O 
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There is a distinction between the two essential singularities we have 
obtained. In this one the singularity occurs everywhere on the line » = 0 
independently of how it was approached, while near Ut x it has an 
essential singularity if Ut + 2+ for fixed ¢ but not for fixed s, no boundary 
conditions being, of course, violated. The reason for this difference in the 
behaviour of the two singularities is in the character of the coefficients of 
6*/e€6n and of 6*ys/elér at the critical values of Ut/x. The first has a simple 
zero when 7» = 0 except that, if in addition € = 0, it has a double zero. The 
second is zero only when ¢ = , having an essential singularity there. If 
the coefficient has a simple zero at the critical value for all y, it is some- 
times, but not always, possible to find a solution of the differential equation 
in question satisfying given boundary conditions. The Rayleigh method, 
discussed in section 3, is an example of the former possibility. 

The results so far may be summarized as follows. After the plate has 
been given its impulsive velocity, the velocity in the boundary layer is 
independent of x until Ut = x. Then, starting at the edge of the boundary 
layer with an essential singularity and moving in towards the plate, the 
x-coordinate begins to affect the solution. Ultimately it becomes the most 
important influence, the change over being about complete when Ut = 2-65x 
and the influence of ¢ dies out like exp(—0-031/%3/23). There may well be 
other solutions of a diffetent kind, but we have been unable to find them. 


7. The effect of compressibility 


The mathematical theory of the boundary layer is complicated by com- 
pressibility effects in two main ways. In the first place, the physical quan- 
tities, viscosity and diffusivity, which were assumed to be constant in the 
incompressible theory, must now vary in a prescribed way with the tem- 
perature. Secondly, four equations, inter-relating the temperature, the 
density, and the components of the fluid velocity along and perpendicular 
to the boundary, have to be integrated, as compared with one only in the 
incompressible case. It is abundantly clear, in fact, that electronic calcu- 
lating machines are necessary to give a complete picture of the flow. 

The literature on steady compressible boundary layers is quite consider- 
able. It may be shown quite easily that when the motion is steady the 
number of dependent variables can be reduced to two, namely the stream 
function and the temperature. To effect further simplification, special 
assumptions must be made as to the nature of the viscosity and the 
diffusivity and also as to the thermal behaviour of the boundary. For 
example, if we assume that the ratio of the diffusivity and the viscosity 
is unity, and that the boundary is thermally insulating, then the tempera- 
ture may be expressed in terms of the stream function, so that the problem 
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e is reduced to the solution of a single partial differential equation of the 
0 same character as that for incompressible flow. If in addition to these two 
u ,ssumptions the viscosity is assumed to be proportional to the absolute 
ry temperature then a suitable transform correlates the compressible boun- 
he dary layer to an incompressible one with a different main stream velocity. 


0 At room temperatures it is known empirically that 


ye iv > T'8/9. (7.1) 
he where 7' is the absolute temperature, and also that the Prandtl number 
pe,,/k 0-715, (7.2) 


where k is the thermal conductivity, and c,, the specific heat at constant 
pressure. The mathematical theory of the steady compressible boundary 

| layer mentioned above can therefore be expected to be qualitatively correct 
ind to exhibit the main features of the flow. The agreement with the more 
exact theory in certain problems is surprisingly good. 


On the other hand, the theory of the unsteady compressible boundary 


: ver is by no means so well developed, chiefly on account of the com- 
th ’ : i ; . : . 
plexity of the equation of continuity consequent on the introduction of the 
+ time derivative. It will be shown here that a suitable transformation can 
in) ‘ : : . ° ° . 
considerably simplify the equations of motion, reducing them to two simul- 
taneous partial differential equations, for the stream function and the 
ibsolute temperature. In the particular problem in which we are interested 
} there is no pressure gradient in the main stream and a complete correlation 
| exists between the compressible and incompressible boundary layers. 
rhe equations of motion appropriate to a boundary layer in the neigh- 
jal 
§ bourhood of a plane wall are as follows 
Tne 
a the equation of continuity, 
en P 
+] cp ( c : : a, 
ne pu) (pv) QO: (7.5) 
cot Ca CY 
ulal ‘ 
the b) the Prandtl dynamical equations 
" ou ( cov Cp C j cou a, 
leu ; p| u 7 I + (1 ) (7.4) 
\ot C2 cy ox oy cy 
sy op ~~ 
ds 0 re; (7.5) 
the Cy 
eal und (c) the equation ot enerey 
ecla [eT oT oT’) cp cp C oT eu\? — 
pe, | u | — ae | Lu . (7.6) 
the | ol Cn cy ot Cx cy\ cy cy 
Fot where x,y are the coordinates of any point measured along and perpen- 


psity | dicular to the wall, and u,v are the corresponding components of the fluid 
pera elocity. In addition there is an equation of state 


blen p = RpT, (7.7) 
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where F# is a constant, and the boundary conditions are 

u v 0 when y 0: u>U,, T+T, ayo, (7.8) 
with a condition on 7' at the wall depending on its thermal character, 
These five equations and five boundary conditions are sufficient to specif; 
the flow. If the flow were steady, @p/ét would vanish, and a stream function 
can be defined.at once from the equation of continuity. For unsteady flows, 


however, we have to use two functions ¢s, Y and specify p, pw, pv in terms 


Cus ; oY = 
pu Po éy > Pp Po éy ’ (7.9) 
Jizt < x,t 


Cus oY - 
py = —Pai- — POs . (7.10) 
OX] 4 ct my 


The relation between p and Y is equivalent to introducing a new coordinate 


of them. If we set 


then it follows that 


Y, instead of y, defined by . 

y | p dy/ po. (7.11) 
0 

The suffix 0 denotes some standard state depending on the problem under 

consideration, and the suffix 1 denotes conditions in the main stream 

(yoo). In terms of x, Y, t 


U +) . (7.12) 
OY} +4 


Cus Cus Cus oY — 
and - = |- —i— ; 5 (7.13) 
Ox yt Ox Yi c } xrt\€ Hy yt 
so that te — (' *) = ( ) ( ") = (" . (7.14) 
Po Cx) ¥4 ¢ } art\OX } yt ot iene 


Hence we may show, after a little reduction, that 


ou. Ou cu ous Cbs eb us 07 ys — 
+4 | 4 SE meee abe Ae Ae (7.15) 
ot Ox =o OYot  eY dxoy ox oY? 
. ou Lp ou 
Further, fh - a 


oy py OY 


so that (1. = = Df - i cs , (7.16) 
gy Y Popo OY\ T eY* 


l cp T 1 op — 
and moreover Zi _ Pi (7.17 
pox Tip, i 


since from equation (7.5) p is independent of y. Hence the equation 0! 


motion is 
“2 a2 2 . ‘4al ‘ . ‘4al a2, |, \ 
Ce ob eb Aub Cys 1 1 CP; , MoPi § (i é "ap (7.18 
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and has been reduced to a similar form to that for the corresponding steady 
8 fow. It is interesting to observe that the significant transformation, 
equation (7.11), is the same as in steady compressible flow, and that the 
important viscosity parameter is not » but »/7', so that those fluids in 
which u is proportional to the absolute temperature are easier to handle 


than those in which it is constant. To complete the transformation, the 








. } energy equation must also be reduced, and in a similar way we get 
] ‘ed c wiv ks ous CT | r CPy cys Cp, 
P el et oY oa ea eY | a | ator &: 
J VoPy ¢ = oT | el (7 19) 
Do OY \p, To « r) Po Po L'\eY?) - - 
10 The equations of motion have thus been reduced to a pair of simul- 
taneous equations in % and 7’, analogously to those for a steady boundary 
layer. Without further assumptions it does not appear possible to effect 
my further real simplification. This is not surprising, since if the fluid is 
_ compressible the unsteady motion in the main stream where the viscosity 
ilies nay be neglected is imperfectly understood. For the problem in which we 
ae | are interested, the pressure in the main stream is constant and the main 
. features of the flow may be exhibited, if we make use of the following 
issumptions which have already proved useful in steady problems: 
- a) the viscosity varies as the absolute temperature, 
b) o | 
jlo) | c) the boundary is thermally insulating. 
| The dynamical equation of motion is now the same whether the fluid 
7.14 s compressible or not and hence the two solutions of the problem are 
formally the same. If the fluid is compressible, however, the solution is 
nly completed when we have transformed back from Y to y. For this 
715 urpose we need the energy integral given by Crocco (5), which is 
T=T)1 on =» U3 u?). (7.20) 
This satisfies the boundary conditions as Y + 00 and as Y > 0, and hence 
(7.16 a 
y= | (14%—* 2(U2-w2y] a. (7.21) 
J | “y F1 
(7.17 ‘ 
Thus one effect of compressibility is to thicken the boundary layer, 
ion cause Y may be taken to represent the incompressible flow and y the 
compressible flow. At corresponding states of the fluids y Y. fUt< 2 
(7.18) | ih, (7.22) 


2\/(vot) 
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and hence when y is large 

y = Y+(y—1)M?2Y(2r9t), (7.23) 
writing M for the Mach number of the main stream velocity. Moreover for 
large F,, yb = UY —2U\(vpt/n), (7.24) 


and hence the outflow v into the main stream is 


v — ue | (%) (7.25) 


If this velocity normal to the plate is applied as a boundary condition 
in the main stream where viscosity may be neglected, a compression wave 
with infinitely sharp front would be found to be radiated outwards. Since 
this wave is inversely proportional to the square root of the Reynolds 
number, its effect on the boundary layer would be small were it not for 
the singularity when ¢ = 0. For an investigation of the interaction of 
boundary layer and compression wave we must consider the full equations 
of motion, which is beyond the scope of this paper. 

When Ut > x the full solution is not known, but the linearized equation 


(section 3) then gives an outflow of 


» = [14% aren en), /(¥ "), ial 
~ TX 


implying a second compression wave when Ut = x. The discussion of the 


boundary-layer equations suggests that the outflow is not discontinuous at 
Ut = x, but that it gradually rises up to its ultimate value as Ut/xr>« 


The author would like to express his thanks to Professor Howarth for 


stimulating discussions on this paper. 
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for THE SOLUTION OF DIFFERENTIAL EQUATIONS IN 
94 THREE DIMENSIONS 


| I. BOUNDARY VALUE POTENTIAL PROBLEMS 





5 By D. N. DEG. ALLEN and 8. C. R. DENNIS 
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a | 
ne SUMMARY * 
ids | } xtension of the use of tl relaxation method from two to three dimensions 
ry is for a long time been desirabl the use of a three-dimensional network has been 
uught to be impractical because it has been supposed that it must require the 
perposition of a set of plane nets each drawn on a separate sheet. An alternative 
ons | = method was proposed by Tranter (1) in which he attempted to replace one three- 
nsional problem by a set of two-dimensional problems; his method is satis- 
tio} factory, but in an extremely limited capacity in that it only solves Laplace’s (or 
Poisson’s juation inside a volume which is bounded by a cylinder (of any cross- 
tion) with plane ends normal to the generators. In this paper we show that 
oO | ided the work is properly set out, then there is no particular difficulty in the 
om se of a three-dimensional network; as illustrative examples solutions have been 
; d for the steady temperature-distribution problem used by Tranter in describing 
L the s method (for « nparison); and also for the electric potential-distribution inside 
US al ladrant trometer. 
» fot 1. THE use of relaxation methods to solve a Poisson-type equation in two 
| dimensions is now familiar? and does not need to be described again. 
Some of the statements and equations required in this description of the 
extension of the methods into three dimensions will be recognized as being 
the counterparts of similar results already proved in two dimensions; it is 


not necessary to prove them again or to elaborate on their significance, 
We therefore only describe in detail here those aspects of the method 
which are either peculiar to three-dimensional problems or have not pre- 


viously been specifically noted in two dimensions. 


2. We propose to show how a numerical solution in three dimensions 

of the equation ) ) 2 
C-UW “W) CO-W . 

. " W (x, y, z) (1) 

EX cy* Cz" 

can be found by determining, to any required order of accuracy, values 

of the required quantity, w, at the joints of a uniform cubical lattice which 

extends throughout the volume of integration; W is supposed to be a 

\ detailed description has been given by Southwell (2). 


{Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 2 (1951)] 
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function of x, y, and z known in so far as a numerical value can be assigned 
to it at all the joints of the lattice. The method consists essentially of 
replacing the differential equation (1) by a set of finite-difference approxi 
mations which relate the value of w at a typical joint of the lattice with 
the values at the six surrounding joints. Fig. 1 represents such a typical 


joint, 0, and the six adjacent joints are denoted by 1, 2, 3, 4, 5, and 6 as 





2 
] 
Pa 

3 o| -* ms 
—— i = 

| 

|s 

Fic. | 


shown; when the mesh-distance between neighbouring joints is taken to 
be h, the finite-difference replacement for the equation (1), holding at the 
joint 0, is known to be 


. Bi tanoss , , ; ee 27 9 
W,+w.+wst+w,+u;+w_,— bw, hv, 0, (2) 


where suffixes denote values of w and W at correspondingly numbered 
joints. In the process of a relaxation computation the values of w at the 
joints change gradually from some initial (and arbitrarily chosen) values 
and eventually converge to the finally accepted values which constitute 
the solution. When this convergence is completed the equation (2) must 








O— 2 —O 


Nic. 2. 


be satisfied at each joint of the lattice as nearly as that is possible, depend 
ing on the number of significant figures which are used; at any stage in the 
process of the computation, before the accepted solution has been reached 
the equation (2) will not necessarily be satisfied at any joint by the w-values 
obtaining at that stage. The amounts by which that equation is not then 
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satisfied are measured at each node by a residual defined by the expression 


F uw u uw “uw 


0 1 2 3 1 Ws pu 


5 —6w,—h? WM. (3) 


The aim of the relaxation method is eventually to make every residual as 
near to zero as possible by continual adjustment of the w-values. The work 
is systematically and easily accomplished by the use of a relaxation pattern 
siven in Fig. 2, which shows the changes occasioned in the values of the 
residuals when a unit increment is added to w at a typical joint of the 


attice. 


3. The major difficulty which has hitherto prevented three-dimensional 
relaxation is that involved in finding a satisfactory technique of recording 
nm paper a calculation which in reality should take place on a three- 
dimensional lattice. Tranter avoids the difficulty by, in effect, using 
Fourier expansions in one (say z) of the coordinates. His solution w is then 
determined in the infinite series form, 

» 4 
w <a @(n)sin nz, 


— 


n=1 
where @(n) is a function of x and y only which can be found by a two- 
limensional relaxation. He thus reduces the one three-dimensional prob- 
lem to a number of two-dimensional problems—to find @(1), @(2), @(3), 
ete.; but this can be done only at the expense of restricting the shape of 
the volume of integration to be cylindrical and by requiring the boundary 
condition on the plane (z constant) ends to be either given values of w 


or of the normal oradient of w 


4. To illustrate how a three-dimensional lattice may be used we take 
Tranter’s example—to determine the steady temperature distribution 
nside a uniform cube when one face is maintained at a constant tempera- 
ture (say 1,000) and all the other faces are kept at zero temperature. The 
equation governing the temperature, w, throughout the cube is Laplace’s, 
le. equation (1) with W 0, so that the formula (3) which defines 


residuals at each joint of the lattice reduces simply to 


F, wy Ws U's U's + Ws + We OW. (4) 


| It may be noticed that we do not need to specify any particular absolute 
size of the cube. If its side is of length L, then its six faces may be taken 
to be the planes x, y, 2 = 0 or L and we will suppose that the face y = L 
is the one on which w is kept at 1,000. The solution was obtained initially 
on a coarse lattice with each edge of the cube containing four mesh-lengths 


i = L/4). This complete lattice covering the whole cube is shown in 
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Fig. 3: boundary values of w are specified at the joints of the lattice jn 
the faces of the cube, and twenty-seven internal joints can be counted: 
from each of these joints six lattice-lines radiate in the axial directions, 
and once initial (guessed) values of w have been assigned at these joints, 
residuals can be calculated by substitution in the formula (4). They may 
then be relaxed by use of the usual methods of relaxation and, in particular, 
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of the pattern of Fig. 2, which has been drawn to correspond with Fig. 3. 
Fig. 3 itself is given primarily to illustrate how a (typical) three-dimensional 
lattice should be drawn in order that a relaxation computation may be 
carried out on it; it can be seen that it allows at each joint of two columns 
of values being written down, one consisting of increments added to the 
w-value there and the other of successive values of the residual. It is 
necessary in drawing the lattice to allow only one set of plane sections to 
overlap one another. To achieve this and at the same time to compress 
the diagram into as small a space as possible necessitates considerable 
relative distortion of the mesh-lengths in the three axial directions. /t is 
this ‘isometric projection’ of the three-dimensional lattice which enables th 
work to be carried out on a single sheet of paper. 

The three dimensions of the region shown in Fig. 3 are apparent in spite 
of the distortion already mentioned; this distortion is in no way trouble 
some because in the actual computations of the relaxation process made 
on it, the lattice is never thought of in three dimensions but instead as a two- 
dimensional network with nodes at all points marked with dots; all other 
(unmarked) points where strings cross are not regarded as nodes of the net. 
This view of the actual working diagram greatly facilitates the relaxation: 


nodes and its effects on the residuals immediately recorded by the usual 
and simple two-dimensional technique of ‘counting strings’. 


in particular a block-displacement can be applied for any chosen group 0! 
| Py ; 
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5. In this example there is obvious symmetry in both of the planes 
y= L/2 and L/2, so only one-quarter of the whole cube need actually 
be used in the computation and in the presentation of results; these are 
given in Figs. 4 and 5. Fig. 4 records values of w, at the joints, as found 
on the coarse lattice (/ L/4). Fig. 5 similarly records values of w found 
at the joints of a finer lattice (/ L8). Although, of course, the actual 





Fic. 4. 


process of relaxation on a three-dimensional lattice is rather more compli- 
cated than it is on a two-dimensional net, it is nevertheless no more 
laborious to arrive at a solution and indeed, with a little experience, 
convergence is probably more rapid in three dimensions than in two for a 
comparable number of joints and nodes respectively. The reason for this is 
that the ratio of boundary to internal joints in three dimensions is greater 
than the corresponding ratio of boundary to internal nodes in two dimen- 
sions; consequently it is easier in three dimensions to liquidate initial 
residuals by the relaxation process of moving them to the boundary, 
where, in a boundary value problem at least, they are not required to 
vanish. Thus in the finer lattice of the example given in this section 
there are 243 internal joints at any of which there may initially be non- 
zero residuals to be removed; but none of these joints are more than four 
mesh-distances from some part of the boundary. A corresponding net in 


two dimensions would cover a square with sixteen mesh-lengths in each 
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side (this net would include 225 internal nodes); and the individual nodes 
of such a square net are anything up to eight mesh-lengths from the nearest 
part of the boundary. 

















6. An exact analytical solution in the form of a double infinite series 
has been found for this problem by Carslaw (3) and is quoted by Tranter. 
Calculation from it of the values of w at the joints of the lattice used in 
the relaxation solution enable a test for accuracy to be made in this first 
example to be solved by three-dimensional relaxation. It is found that for 
the coarse lattice of Fig. 4 the maximum error in the w-values is 22, 01 
2-2 per cent. of the greatest w-value. For the finer lattice of Fig. 5 this 
maximum error is reduced to 1-7 per cent. No particularly useful purpos 
would be served by making any comparison with Tranter’s results, and 
this is not done since it would entail considerable extra computation: fot 
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he gives solutions for the functions (1), ete. (cf. section 3), but uses them 
to find w itself at only one} point. 


7. In general the bounding surface of the volume of integration may be 
curved and then some of the strings of the lattice radiating from joints just 


inside the boundary would be cut short by it. Fig. 6 illustrates a possible 





2 
5 
” 3 = ee 


case where the string OB actually meets the boundary at the point B. 
The point 1 then would lie just outside the volume of integration and any 
value w, which may be attached to it is a fictitious value; the residual at 
the point 0 is still given by the formula (4), but in it is contained this 
fictitious w, which must be eliminated before the residual Fj can be 
evaluated. The necessary elimination can be performed by making use of 
the fact that the boundary condition provides a known value wz. If OB 
is a fraction € of the complete mesh-length h, then it is easily shown that, 


with an error of O(h?), 


2 1—é 2(1—&) (5) 
Th uw w - -Ww » 
1 B aes 0° 
€(] 1+€ 3 


For with 0 as origin w may be expanded in powers of x along the line 301: 


uw Wo per ; qu? Tr eeey 

where p = (€w/éx), and q €*w/éa*),/2!. Putting 2 equal in turn to €h, 
h, —h in this series, we get 

WR Wot péh I gé*h*, 

Vy Wo + ph i qh*, 

Us wW,y— ph+-qh*, 
with neglect of terms containing h*® and higher powers of h. Elimination 
ol p and q from these three equations yields the relation (5). When this 


le. at the centre of the cube ; it may be remarked that the exact value there is obviously 
1000/6 since the effect of maintaining one face of the cube at a temperature of 1,000 is by 
mmetry one-sixth of the effect of so maintaining all six faces. Any numerical method 
ch makes use of this symmetry must yield a correct value at that point; and testing 

h a solution there affords no check of accuracy whatsoever. 
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relation is used to substitute for the fictitious w, in the residual formula 
(4) we find that it becomest 


» 


») 9 
F 1b. 1-68 ogee, A. = ap a W rp At —lw 6 
0 2 1 i ‘ is") se B ) 0° (0) 
; l g &( l T £) > 


A short string radiating from a joint 0 in any of the other five directions 
can be dealt with in an exactly similar manner; and an obvious combina- 
tion of two or more of these cases can be made when more than one of 


the strings radiating from a joint are cut short by the boundary. 


8. A modified residual formula is also used at joints which lie in a plane 
of symmetry. Thus if 0 typifies such a joint and (Fig. 1) the plane of 
symmetry is the plane 012, it is not necessary to record the field on both 
sides of this plane: the joint 6 (say) would not be shown since it is simply 
the reflection of the joint 5, consequently in the residual formula (4) the 
value w, is in effect fictitious; it is eliminated by means of the relation 
We, (7) 
so that the residual, F,, is given by 

Fi WwW, +W.+Ws+uU,+ 2w;— bw. (8) 
It is worth mentioning that one can satisfy a normal gradient boundary- 
condition at a plane boundary by this same treatment; for a plane of 
symmetry is in effect a boundary (of the volume which is recorded) with 
a condition imposed there that cw/év = 0 (v being the direction normal to 


the plane).{ If, more generally, we require to impose the condition 


éw/ev = N(x, y,z), then the only difference is that (7) is replaced by 
’ : ») y 
We w,+2hN,, 

and (8) by Fy, = w,+w.tw,+w,+ 2w;—6wy+2hNp. 


9. As our second worked example we take one which includes a curved 
boundary and therefore a use of the irregular star residual-formula (6), and 
for which also Tranter’s method is not directly applicable. The problem 
is to determine the electric potential distribution inside a quadrant electro- 
meter whose dimensions are shown in Fig. 7. The region of integration is 
bounded externally by a short closed circular cylinder divided axially into 
four quadrants, and internally by a plane vane, parallel to and midway 
between the ends of the cylindrical box. The shape of this vane is shown 

This derivation of the ‘irregular star’ residual formula differs from that given in tw 
dimensions by Southwell (2), because one cannot, in three dimensions, make use of th 
analogy which interprets w as the transverse deflexion of a uniformly tensioned string 


network. 


t It should perhaps be mentioned that across a plane of symmetry (7) is an exact relatiot 
whereas across a true boundary it is only the usual approximation to the condition éw/cv = 9 
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EQUATIONS 


in the plan view (Fig. 7a); its position along the axis of the box is shown 
in the elevations (Fig. 7 b and c). 
The equation to be satisfied by the electric potential V, in the region 


between the vane and the box, is again Laplace’s equation, V7V = 0. The 
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Fic. 7 
undary conditions given are that the potential is maintained at pre- 


scribed values on the vane (J) and on each pair of opposite quadrants of 


the box (V, and V3). In computation these values were taken to be 
on the vane, | 200, | 
on one pair ol quadrants, Vy, 4, (9) 
on the othe pair of quadrants, A J 


because the thickness of the box is small compared with the diameter of 


| lls cross-section, a smaller mesh-length (h/4) was used in the direction 
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parallel to its axis than the length (4) in the two perpendicular directions; 
this circumstance implies a modified residual-formula, 

Fy, = V.th+h+%,4+ 16+ 16%—36h, (10) 
and a correspondingly modified relaxation-pattern. On the finest lattice 
used, there were four of the smaller mesh-lengths in the half-thickness of 
the box and eight of the larger mesh-lengths in the radius. The results 
are shown in Fig. 8 which records values, at the lattice-joints, of the 
potential, V. Conditions of symmetry make it sufficient to give results 
over one-quarter only of the whole volume. 
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SUMMARY 


Relaxation methods have been applied in recent years to obtain solutions in two 


ensions to various second-order linear partial differential equations. An example 
is provided by Laplace’s equationt which, though simple, is typical of all the equa- 
ns which have been solved in that it is of the so-called ‘elliptic’ form and in that 
| solutions which have been found have been required to satisfy a single boundary 
condition at all points of a closed boundary. There are many problems, however, for 
vhich the relaxation method has not hitherto been applicable because either one or 
isually) both of these properties do not hold. Typical examples of these are provided 
when solutions are required of the heat-conduction equation (which is of ‘parabolic’ 
form) and of the wave equation (which is of ‘hyperbolic’ form) ; and also, occasionally, 
f Laplace’s equation itself (which, although of ‘elliptic’ form, may not be associated 
vith a single boundary condition at all points of a complete boundary). 


In one dimension inherently the same difficulty distinguishes what are usually 
snown as ‘jury’ problems from ‘marching’ problems. A method of overcoming these 
fficulties is presented here ; the method was originally devised in order to bring the 
heat-conduction equation within the scope of relaxation and this paper is principally 


concerned with that application; the underlying idea is, however, introduced, for 


implicity only, by illustrations in one dimension. 
1. RELAXATION methods were first devised, and for many years developed, 
as an alternative to classical analysis in order to find particular solutions 
of differential equations. By now, however, their procedures have been so 
lar extended as sometimes to influence the general plan of attack. One 
example, illustrated by Fox and Southwell (2), in the determination of 
extensional stresses in a flat elastic plate, is the more convenient employ- 
ment of two simultaneous differential equations instead of an orthodox 
reduction of them to a single equation of higher (i.e. biharmonic) order. 
That tendency is maintained in this paper by the introduction, solely on 
jrounds of relaxational convenience, of operations which would not be con- 
templated in classical analysis because they would make the problem not 
easier, but in fact harder, to solve. Thus, in order to make the (heat- 
conduction) equation, 


ov on1 
] 
ot Ox? (1) 
\ very full discussion has been given by Southwell (1). 


[Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 2 (1951)] 
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tractable by relaxation, a function w is introduced which is related to » 
by the equation a oe 

_ Ow Bw, 
tT Gaz? 2) 
elimination of v between equations (1) and (2) leads to a modified problem 
governed by the equation aw otw 
a ex! ") 
with, of course, appropriate boundary conditions. Clearly, equation (3) is 
harder than equation (1) to solve by the methods of classical analysis. In 
contrast, the device just outlined makes this problem for the first time 
soluble by relaxation. An example is worked in detail in sections 9-10 with 
results which are presented in Fig. 1. 


2. The difference between ‘jury’ and ‘marching’ problems was discussed 
by Richardson (3) who gave them those names. In one dimension a march- 
ing problem is one in which the end-conditions to be imposed are all stated 
at one end of the range of integration of the differential equation to be 
solved; and numerical solutions are usually constructed by means of a 
‘step-by-step’ integration process. In a jury problem, on the other hand, 
the conditions are imposed in equal numbers at both ends; classical analysis 
then must have recourse to a synthesis of solutions which in numerical 
work are found, successively, by a marching process; as such they are 
liable to accumulated error, and the fact that the relaxation method 
enables jury problems to be solved directly is one of the features which 
gives it claim to importance. (The problem being, so to speak, ‘tethered’ 
at both ends, or—in two dimensions—all round a boundary, the solution 
cannot ‘take charge’ anywhere in the region.) 

[t is not the intention of this paper to suggest that ‘step-by-step’ 
methods of integration are not perfectly satisfactory to solve many kinds 
of marching problems, nor is it suggested that the method described here 
necessarily obtains more accurate results; only experience will in time 


indicate which method is the more suitable to employ for given types of 


problem. In particular the numerical examples given as illustrations in 
this paper certainly yield satisfactorily to a ‘step-by-step’ attack—indeed, 
for that matter, exact solutions can be found for them by analytical means; 
they have been chosen here as intrinsically simple illustrations of an 
alternative method which is being developed to solve much more compli- 
cated problems, some of which present great difficulties to an approach 
either by orthodox analysis or by ‘step-by-step’ operations. 


3. Itis, in fact, only to problems of jury type that the relaxation method 


is directly applicable; and the significance of the procedure outlined in 
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1G. 1, The distribution of temperature (v) in a rod when one end is suddenly main- 
ained at v 


100, the other end being kept at the initial temperature of the whole 
0. (a) The variation of temperature with time (t) at five equally-spaced 


nts of s 


ubdivision along the rod ; (6) the variation of temperature with position (2) 
at instants of time, ¢ [*/144x«V2)(1, 3, 5, 7, 9, 11). 
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section 1 is that it makes a jury problem of one which as presented was 
of marching type. It was devised in consequence of purely relaxational 
considerations of the circumstances which are necessary for solving a 
differential equation. Its effect would be nugatory in a classical attack, 
which would only yield, for instance, as a solution to equation (3): 


w ] [ | [ V dex} da+W, (4) 


where V denotes the general solution of equation (1), and W the general 
solution of equation (2) when v is put equal to zero. Then from (4), accord- 
ing to (2), it follows that » = V, simply, is the solution to equation (1) so 
that the introduction of the function w has achieved nothing. The utility 
of the procedure from a relaxational standpoint is only revealed when the 
differential equations are replaced by their finite-difference approximations. 

4. For clarity it is well to consider firstly equations in one independent 
variable (x) only; a relaxational solution of such an equation consists of 
numerical values of the required quantity or dependent variable (v) ata 
number of points of subdivision equally spaced along the range of integra- 
tion. In terms of v-values at such points, finite-difference approximations 
to low order v-derivatives are well known: Fig. 2 represents five typical 


<_—_—h- 
eae See es 
11 3 O { ‘ a 
Fic. 2. 


successive points of subdivision, and the approximations to the first four 
derivatives at the central point 0 are 
2h(dv/dx)y = V1—Vs, 
h?(d?v/dx*), = v,-+03—2v9, " 
2h3(d3v/dx*), Vyp— 20, +203—1y,, | 
and h4(d4v/da*), —4v,— 4054-699 +?y), 
where suffixes denote values obtaining at correspondingly numbered points, 
and h is the length of the interval between successive points. It may be 
noticed that the approximations to the first and second derivatives contain 
values of v at points at distances up to, but not exceeding, one interval 
from 0; we shall say, briefly, that these approximations ‘extend’ over one 
interval on either side of 0. It is evident that, in general, the approxima 
tions to the derivatives of orders 2n—1 and 2n both similarly extend over 
n intervals on either side ot the central point. 
If, then, the governing equation is of order either 2n—1 or 2n, its finite- 
difference approximation, holding at the point 0, extends over intervals 
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on either side of 0; and thus when 0 represents either end-point of the range 
the finite-difference approximation involves n ‘fictitious’ v-values at n 
fictitious points outside the range. As an essential preliminary to the use 
f relaxation, these fictitious values have to be eliminated; thus  condi- 
tions are required at each end that this elimination may be achieved. In 
, jury problem this requirement is satisfied and a relaxational solution can 
be directly effected. But in a marching problem, as presented, no conditions 
wre available at one end—we have instead some unwanted conditions 
additional to the number that are wanted) at the other end; a relaxational 
technique cannot then be directly applied. 

Since we have just seen that relaxation methods can only be utilized to 
solve a problem when it is expressed in jury form, and that then the order 
if the equation is double the number of conditions available at either end 
f the range, it follows conversely that if m end-conditions are required to 
e imposed at one end, then the governing equation—if not of order 2m as 
presented—must be transformed into a modified equation of that order 
before a relaxational technique can be used. In general m conditions may 
ye given at one end and p conditions at the other, and it is sufficient to 
suppose that p < m; the order of the governing equation as presented 

ould then be m+-p. If p = m the problem is of jury type and needs no 


further discussion. If p < m the governing equation must be transformed 
so that its order is increased by m—p to become 2m; then it will be possible - 


to impose m conditions at the one end while at the other end the original 

conditions must be augmented by m—p extra conditions and the prob- 
em will thus be converted into jury type. As an extreme case, in a 
iarching problem m conditions are given at one end but none at the other; 
he transformation must now double the order of the governing equation 


rom m to 2m. 


5. The transformation of the governing equation may be accomplished 
means of a change of the dependent variable by substitution in the 
juation, in place of v, a differential function of a new dependent variable, 


We make no attempt to lay down rules according to which this substi- 


ossibilities which would suffice and a compromise must be made, firstly 
tomake the substitution itself as simple as possible, and secondly to make 
the transformed equation as convenient as possible for solution by relaxa- 
ion. It is apparent that the order of the differential function of w which 
must be substituted for v is equal to the desired increase in order (i.e. to 

p) of the transformed equation above the original equation; and hence 


[hus keeping to a minimum the work which has to be done after the transformed 
juation has been solved 
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is also equal to the number of extra end-conditions which have to be 
supplied. These extra end-conditions can be arbitrarily chosen provided 
only that they must not imply any restriction on v or v-derivatives: 
normally they would consist of assigned values to w and the first (m—p—1) 
w-derivatives (such conditions can manifestly never be interpreted in 
terms of v). 

The transformation of the governing equation could be accomplished by 
an alternative method which does not require the substitution for v of a 
differential function of a new variable, w, but instead requires us to operate 
through that equation with a differential operator so as to increase its order 
by the desired amount. This differential operator has to be chosen with 
aims similar to those which govern the choice of the substitution in the 
first method, and in the worked examples of this paper there would be 
little to choose between the two methods. We shall, however, adopt the 
first method here, partly because it is perhaps the less obvious of the two 
and therefore the more desirable to emphasize; and partly because it has 
been found to be the more convenient method to use for more advanced 


problems. 


6. As an example to illustrate the procedure we consider the equation 


d*v/du? = f(x), (6) 
to be solved over a range of integration a < x < 5b; f(x) is a function of x 


which is supposed known only in the sense that its numerical value can be 
calculated or measured for any value of x within that range. In order to 


define a particular solution of (6) two end-conditions must be given; if 


these are stated one at each end of the range, an immediate and straight- 
forward application of the relaxation method will yield the solution to any 
required degree of accuracy. But if the conditions are both given at one 
end of the range, as for instance, that at 

z= 4, v= & 


and dvujdx = K, 


then the problem, as presented, is of marching type and the equation must 
be transformed into another of double the original order before relaxation 
can be applied. 
A simple and sufficient substitution is given by the relaxation 
v = d*w/dz?, 
which transforms (6) into the modified equation 


d‘w/da* = f(x). 
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The given end-conditions (7) become, in terms of w, that at 


(10) 


x a, d?w/dx? =k ! 
and Bw/da? — K, 


while the two extra conditions to be supplied may be chosen, as simply as 
possible, to be that at 


x b, w 0 ! at 
) 
and dw/dx 0. 


The problem is now converted to be of jury type, requiring only the 
solution for w of equation (9) subject to the end-conditions (10) and (11); 
it can therefore be solved by relaxation. And when w has been so found, 
>can be at once derived by use of equation (8). It should be mentioned 
(cf. section 3) that any other choice for the extra conditions (11) would 
only alter the w-values by a quantity A+ Bax (where A and B are some 
constants) and would therefore not make any difference in v. 

If the alternative (second) method of transforming the governing equa- 
tion were used for this example, it would require (6) to be replaced by the 


f"(x); 


the end-conditions (7) would be unaltered, and the extra end-conditions 


modified equation d4y/dat 


at « = b would comprise (6) itself, together with 

d3v/dx? = f'(a). 
In this method the extra end-conditions are, of course, no longer arbitrary, 
but must instead be derived from the original governing equation. 


7. As a second example in one dimension we take for solution the first- 


dw 
dx 


order equation 


over the range 0 < x and with the end-condition given to be that at 
r 0. v :. (13) 
\gain the order of the equation must be doubled by a preliminary trans- 
formation. In this case an appropriate substitution is 
dw/da+-w; (14) 
thus (12) is transformed to become 
d*w/dx?—w—-ax 0, (15) 
and, in terms of w, the end-condition (13) becomes that at 
== 2: dw/dx+w b. (16) 
The extra condition is chosen to be simply that at 


s=1, we @. (17) 
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The problem is now converted into jury type—to solve the equation 
(15) for w subject to the conditions (16) and (17). The solution for » 
(multiplied by 1,000 to avoid decimals) is shown in Fig. 3; it was easily 
effected by the normal relaxation procedure in less than 30 minutes, with 
the range divided into ten intervals. Fig. 3 also records (in brackets) values 








| | 
| 
(!000) a ail 1.241 (1.397. (1.579) (1.792)} (2.037) (2.319) (2.641) (3.007 (3.422) 
-36450 ~31987| aH -23939} -20253 -1673 -1332 -9991 -6688 ~3372 ° 
x=0 xml 


Fic. 3. 


of v computed afterwards from the finite-difference approximation to (14), 
Vo O(W1— Ws) + Wp. 


These values of v have been compared with the exact values (derived from 
the analytical solution v = 2e*—x—1); they are nowhere in error by as 
much as } per cent. They are also negligibly different from the values 
which are obtained when the alternative (second) method of transforming 
the governing equation, which is suggested in section 5, is employed. 


8. Having explained and illustrated the methods of this paper in one 
dimension we now proceed to apply them in two dimensions to the solution 
of the heat-conduction equation 


© aie, as) 
ot Ox? 
where v, a function of time, ¢, and distance, 2, measures the temperature 
in a linear flow of heat; « is the diffusivity of the medium. 

Two difficulties present themselves if an attempt is made to solve this 
equation directly by relaxation. A rectangular net would be used to cover 
the region of integration, and Fig. 4 shows a typical group of seven nodal 
points on such a net; the mesh-lengths in the x- and ¢-directions are denoted 
respectively by h and H. The finite-difference replacement of equation 
(18), holding at the typical node, 0, can be written down by use of the 
approximations (5), and leads to the definition of a residual there by the 
er Fy = 1, +03—2vp—(h?2/2«H)(v.—14). (19) 
If we conveniently choose mesh-lengths so that h? = 2«H (as we clearly 
may, cf. section 9), then the residual formula (19) yields a corresponding 
relaxation-pattern (illustrated in Fig. 5), which shows the changes made 
to the residuals when a unit increment is added to the value of v at 4 
typical node of the net. 








Th 
easily 
incre! 
total o 


Thus: 
to anc 
moun 
(ineid 

The 
liming 
equat 
of the 
t-deriy 
the tr 


























NON-ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS 217 
yn The first difficulty is apparent in this relaxation pattern; residuals are 
a easily moved about in the x-direction, but it may be observed that an 
ly increment added to a v-value at any node makes no change at all to the 
| total of the residuals on the mesh-line through that node in the x-direction. 
es , 

* 
4 f2 
| H 
+ h-—-— 
—e@ Y * ad ene. 

11 3 e) 1 9 

y ®4 





ca Fic. ov. 


| Thus the pattern does not easily allow residuals to be moved from one node 
the to another in the ¢-direction. This difficulty would not, in itself, be unsur- 
the mountable; we take notice of it here only because we shall have later an 
incidental) opportunity of removing it. 

The second difficulty cannot be overcome, however, except by a pre- 
arly | liminary transformation of the governing equation (18); it is that the 
ling } equation is of marching type in the time coordinate. The principal object 
ad of the transformation is therefore to double the order of (18) in the 
at a -derivatives, but at the same time it is chosen so far as possible to make 


the transformed equation both simple and convenient for the application 
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of relaxation; the substitution 


: ow, Ow 
ES Gat? (20) 


se on @, (21) 


and the finite-difference approximation to this equation, on a rectangular 
net of sides h and H in the z- and f-directions, is 


(h4/x?H?)(w.+w,—2w,) +4w, +4w,—6w,—w,—w,, = 0. (22) 


Fic. 6. 


The relaxation pattern found from this approximation does not suffer from 
the disadvantage which was apparent in that of Fig. 5; we can indeed by a 
suitable choice of the mesh-interval, H, make the pattern derived from 
equation (22) such that the effect of an increment added to w at a node 
of the net is to move the residual away from that node equally in all four 
directions. To do this requires that we choose 
H = h?/ew2, (23) 
and then equation (22) leads to residuals being defined by the expression 
Fy = 2w.+2w,+4w,+4w,—10w,—w,— wv). (24) 


The corresponding relaxation pattern is shown in Fig. 6. 


9. We take as an illustrative example the determination of the tempera- 
ture, in a medium which extends from x = 0 to x = L, which is initially 
at zero temperature throughout; where the end x = 0 is kept at zero 
temperature, but where at time ¢ = 0 the other end 2 = L is raised to and 
subsequently maintained at a temperature of 100. In terms of w this 
problem requires the solution of the equation (21) over a rectangular region 
in the (2, t)-plane bounded by the lines x = 0, x = L, t = 0, andt=T; 
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T measures the extent in time for which the solution is required. The 
boundary conditions to be imposed are that 

() ‘ 

ow Ow 
on x Q, —+t+k— = 9, 
ac ox 
ow C2w | 
) on x L, —+«— = 100, 4 (25) 
ot Ox" 
, ow Cw 

” and ont = 0, —+ce_ = 0, | 

ot ox* 

2 When a solution for w has been found satisfying equation (21) and the 
conditions (25), then v, when calculated from the relation (20), must 
automatically satisfy all requirements. But to find a solution for w we 
need some extra boundary-conditions—one on each of the lines x = 0, 

| ¢g=L,andt T’. These are arbitrarily chosen to be: 

on 2 0, w 0, 
on 2 L, wo = 0, (26) 
and on ¢ T, dOw/ot = 0. 

In choosing these extra conditions we have only to ensure that they cannot 
be interpreted as implying any condition on v or its derivatives, and that 
they are everywhere consistent with one another and with the other con- 
ditions (25). It follows, incidentally, that, by virtue of their being so 

- chosen, the conditions (25) therefore simplify to become: 

Va a ai: as 

on x VU, O-W/ Ox" 0, 

yn 

de on x L, K 6?w/dx* 100, (27) 

Y | and on t= 0, Cw/ét+Ké?w/dx? = 0. 

10. In finite-difference terms the second of the conditions (27) becomes 

99 

rary) 

| (« h*)(w, W,— 2Wy) 100: 

10n 
so that it is convenient to work numerically not in terms of w but in terms 

24 , , om ; ; 

“| of a multiple cw/h? (= w’, say). Then the relation (24) is replaced by 

Fy = 2w.+ 2w4+ 4u}4+ 4u3—10wy—w,— wy, (28) 

Ta and the finite-difference replacements for the boundary conditions (26) and 

ully (27) are that 

| 

ero ; ’ , 

on x = 0, wy = 0, and w,+w, = 0, 

wna 

hin on 4 [, Uy = 90, and w)+w;, = 100, 

nl (29) 

9 


: @, 


+ (w+ w— 2) 
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These conditions are just sufficient for the elimination of all the fictitious 
values which occur in the use of (28) in order to calculate the residuals. 
We can thus find w”’ as a purely numerical quantity and, when it has been 
evaluated, v can also be found, without specifying any particular values 
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for L and x; for, replacing (20) by its finite-difference approximation, using 
(23) and substituting for w in terms of w’, we find that 


, 1 , ¢ J Jf J le 
9 = Wy +we— 29+ (we—U4)/V2. (30) 


The solution is, of course, not independent of the values of L and x; but 
these values only fix the time-interval of the rectangular net at the nodes 
of which v is determined. This interval is given by equation (23); when 
n mesh-lengths are used in the x-direction it follows from that equation 
that the mesh-length in the time-direction is L?/n2«V2, and therefore that, 
if NV mesh-lengths are used in this direction, the total time-range covered 
by the solution is 7 = N L2/n®«v2. 
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The solution has been obtained on two nets of different sizes. In the 
first, or coarser, net the space-range, J, was divided into six mesh-lengths 
(n = 6) and the time-range extended over nine mesh-lengths (N = 9). 
Thus the total time-range was L*/4«v2. It has not been thought necessary 


to record values of w’ and Fig. 7 presents simply the nodal-values of 
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deduced by use of equation (30). These results are presented graphically 
in Fig. 1 in which the curves show, for different values of x, the variation 
of temperature with time. 

On the finer net twelve mesh-lengths were taken in the x-direction and 
also in the ¢-direction (n N 12). In this case the total time-range was 
only L*/12«V2, but within that range much greater detail of the solution 
is obtained. The results are recorded in Fig. 8. From these values curves 
showing the variation of temperature with x, for different values of time, 
were drawn and are also shown in Fig. 1. 

This particular example of a solution of the heat-conduction equation 
can, of course, be worked analytically by a use of Fourier series; it has been 
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chosen here only to illustrate the methods of this paper. The exact solution 
allows, however, of a check of the accuracy of the numerical method: the 
smoothness of the curves obtained and shown in Fig. 1 indicates that it is 
unnecessary to check all of the values given in Figs. 7 and 8; and in a 


number that have been checked the errors do not exceed 3 per cent. of 


the maximum temperature. 
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TABLES OF TWO INTEGRALS AND OF SPIELREIN’S 
is INDUCTANCE FUNCTION 

a By ALAN FLETCHER 

(Department of Applied Mathematics, The University, Liverpool) 


[Received 24 October 1950] 
SUMMARY 
Tables are given for a 0(-01)1 of the integrals 

1 1 

| (K—E)dk and | (K—E)k- dk, 

x x 
where K and E are the complete elliptic integrals with modulus k, and also of an 
expression due to Spielrein which involves both of the above integrals and occurs in 


nnexion with the self-inductance of thin disk coils. Tables are also given of certain 
auxiliary functions and of coefficients required in their calculation. 
1. Account of previous work 
TuIs paper is mainly concerned with the use and extension of formulae of 
known type; no adequate connected account being available, a concise 
summary is necessary. 

The functions to be tabulated arise in connexion with the self-inductance 
of thin disk coils. If the windings consist of a very large number, J, of 
umost circular turns evenly spaced in one plane between two concentric 


limiting circles of radii a and A, the self-inductance in absolute electro- 








magnetic units may be written as L = N?A f(a), where f(a) is a function 
of the single parameter « = a/A 1. A ‘current-sheet’ assumption of the 
usual kind is made; correction for thickness of wire, finite spacing between 
| turns, and so on, is possible but is neglected here. With such assumptions, 
\) is a definite function, but its mathematical determination encountered 
difficulties, and only approximations were at first available. An exact 
expression was, however, published in 1915 by Spielrein (1) and also in 
1919 by Butterworth (2); this expression, which will be called Spielrein’s 
inductance function, may be written, introducing the notations J and J 
lor convenience 16a 
fla (I—o3J), (1) 
3(1— a)? 
: ; 
where | | (K—E) dk J | (K—E)k-> dk, (2) 
x 
K and E being the complete elliptic integrals of the first and second kinds 
with modulus k. From his exact formula Spielrein deduced series approxi- 
mations both for small « (say a 0-5) and for large « (say « > 0-5). He 
also constructed some useful but sketchy tables of f(a). 


[Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 2 (1951)] 
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For the larger values of « an excellent approximation, found by a differ. 
ent method, had already been published in 1914 by Lyle (3), extending an 
earlier approximation due to Rayleigh and Niven; see (4) and (5). Correct- 
ing a known error (p. 435, for 4298579 read 98579), Lyle’s expression, if 
rewritten as a formula for f(«), with Lyle’s c/a replaced by 2y, becomes 


' oe 
f(a) 2m(1-+-a)(SIn So}, 3) 
z 
where 
wher ' ee 103. 
S, 2 2. y?+- -»/4 .o 98, 
24 2880 " 98 3.6.7 
(4 
4! 985 
g —~1_ 42 1 wore 
> | O96” 150° 213 3252_72 
l—oa radial width . 
and y z (5) 


ita mean diameter’ 
This appears to be the best previously published approximation for large a. 
[t should be added that Lyle’s paper was of wider electrical scope than 
Spielrein’s; he investigated coils of mean radius a and rectangular section, 
the axial dimension of the rectangle being b and the radial dimension ¢ 
Thus formulae (3, 4, 5) relate merely to a limiting case, b 0, in his work, 

Formulae similar to (3, 4, 5) were deduced by Spielrein from his exact 
expression by an application of Landen’s transformation. He carried the 
approximation to a lower order than Lyle, in that he did not obtain the 
terms in y® in (4). On the other hand, it was Spielrein’s exact expression 
which enabled the limiting case of a disk coil to be seen in its proper mathe- 
matical perspective, and his work is basic for the present paper, which adds 
terms in y® and y". 

Some details of the effect of the Landen-type transformation (5) will be 
required. In the J, J notation used here, Spielrein’s results, which are 
easily checked by the usual Landen formulae for complete integrals (6), 
show that 

z o\ 17 . . 
pa { ty )A—28 dy, J Lit 2 allie (6) 
(1 | y)8 J (i y)3 
0 0 
where K and FE now have modulus ,/(1—y?). On using the usual expressions 
(7) of the form 


K = K,In(4/y)—K,, E = E,\n(4/y)+ £,, (7) 


where K,, K,, E,, E, are power series in y?, and performing on series the 


various operations demanded by (6), one obtains formulae, essentially due 
to Spielrein, which will be written in the form 


I = L,ln(4/y)—4,, J = J,\n(4/y)—A, (5) 
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where 
3 25 43 1057 . 1507 
] as aoe i oe , 64, 
>. 12° 16’ ' 320” — 384” , 
(9) 
.. 2. Bs 16131 . 6835 ,. 
a a v4 aye at -y° — “Tm ce0g 
. ‘ 4 Q ° 64° ~— 3200 1152 
and J,, J, have the same terms with all signs positive, i.e. 
Jy) I(—y), J,(y) —I,(—y). (10) 
Terms involving y’ to y!* are added to (9) below. 
Using (5), equation (1) may be written 
(1+y)l—(1l—y 8 J 
T(x) 2r7r(1 x) — Y) C-yY) —~. (11) 


oy 


Substituting (8, 9, 10) in (11), one finds, with Spielrein, Lyle’s formulae 


3, 4) with the y® terms of S, and S, omitted. 


Integrals of comple te elliptic integrals 

In describing previous work, it remains only to consider the type of 
integral occurring in (2). The integrals J and J could be replaced by others, 
using relations mentioned by Spielrein and examined in detail by Miiller 
8); see also Jahnke and Emde (9). In virtue of the known differential 


relations dK E K dE : E K 


5 es 12 
dk k(1—k?) k dk k — 


it is shown by Miiller that recurrence formulae exist which, in general, 
connect | Kk' dk with [ Kktdk, and [{ Ek dk with ( Ekt-? dk; but 
these recurrence formulae give no connexion (owing to vanishing of 
coefficients) between { K dk and | Kk-? dk, nor between [ Kk dk and 
Kk-1 dk, and similarly when K is replaced by £. Moreover, 


2(Edk=|Kdk+kE,  [ Ek+dk=[ Kkodk+#. 


Thus, restricting ourselves to integral indices, the results are that | Kk" dk 
ind | Ek’ dk may be expressed rationally in terms of: 


(l) kK, Hifr 1, 3, 5,... or —2, —4, —6.,..., 
(2) k, K, E, | Kdkifr = 0, 2, 4...., 
(3) k, K, BE, { Kk-dk ify ; wh, Be. 


Hence all integrals of the form { Kk" dk and { Ekt dk, where r is integral, 
may be evaluated in terms of, say, ( K dk and ( Kk-' dk, together with 
k, K, E. 

14 Q 
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. ; 
We also have | K dk = 24, | Edk = G+4, (13) 
0 0 


where G is Catalan’s constant, defined by 





47 1 
i; : 
G = $e = | tan 4H 
2) sind , t 
0 0 
i—3 245 27-21... = 0-91596 55941 77219 (14) 
It follows from (2) that 
1 
1(0) | (K—E) dk = G—4, (15) 
0 
and from (1) that 
f(0) = §a(2G—1) = 6-96957 04257. (16) 


In considering the whole class of integrals of the above kind, it might be 
simplest to consider { K dk and { Kk-1dk as the basic integrals in terms 
of which the others may be expressed. But the pair J and J which occur 
directly in Spielrein’s formula have themselves an elegant mutuality visible 
in equations (10), and will be tabulated below. In terms of them, 
we have 1 


| K dk = 21~—a8+1, 





. 


; 

| E dk I—a«E+1, 
_ | (17) 
| Kk-1 dk 2J—(1—a?)a"*K+aH—1, | 


x 





1 
| Hk- dk = 2J—(1—o?)a-*(K—E), | 
x J 
where the complete integrals K and E on the right have modulus «. Thus 
it would be easy to tabulate other integrals of the class if desired. 


2. Remarks on previous work 

The excellent pioneer analysis of Spielrein was accompanied by two less 
attractive numerical tables of f(x). Firstly, he gave f(a) to 4-6 decimals 
for « = 0(-05)-9. Only 4 decimals are given in the middle of the range, 
where values are hardest to calculate. The sixth decimal, when given, is 
always wrong by from 3 to 8 units, and is thus almost worthless; even the 
value of f(0) is wrong by 3 units of the sixth place, although Spielrein gives 
G to 11 decimals and thus has material to determine f(0), by (16), to about 
10 decimals. Secondly, he gave f(x) to 1-3 decimals for « = -05(-01)-95, 99. 
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The first decimal is wrong at « ‘99. It will be noticed that arguments 
are omitted, presumably as unimportant in practice, near 0 and 1. The 
erratic number of decimals given makes interpolation troublesome. 

The case of a thin disk coil is admittedly less important in electrical 
practice than the well-known and mathematically easier case of a thin 
cylindrical coil. Nevertheless, it has its place; thus it is noticeable that 
Lyle (3), wishing to check his formulae for coils of rectangular section, was 
able to do so only at the cylindrical limit (¢ = 0), and not at the disk limit 
(b = 0), because no expressions for disk coils better than his own were then 
available. (The present paper will confirm that Lyle’s coefficients as 
amended in (4) are correct as far as they go.) Spielrein’s table only par- 
tially fills the gap. Grover in his recent book (10) has tabulated 

P = 2f(a)/(1+«) 
with argument y; this variant form has considerable merits, and the table 
is useful when 5 figures suffice. For details of this and a few small tables, 
see Fletcher (11). A recomputation of f(«) seems not out of place. 

There is a further, and very strong, reason for recomputation. As f(a) 
depends on two integrals, it is impossible to convert any table relating to 
disk coils so as to obtain numerical information about either integral 
separately. It seems a pity to refrain from tabulating two expressions 
which would bring a whole new set of integrals into the class of functions 
which may be regarded as known for computational purposes. The inte- 
grals J and J, or some other basic pair, should be tabulated, as well as an 
inductance function. 


3. Aim of present paper 

The present paper aims at tabulating J, J, f(«) for « = 0(-01)1. A reason 
additional to those given above stimulated the author to undertake the 
necessary labour. In 1940 he constructed (6) a 10-decimal table of K and Z 
for k 0(-01)1; this arose out of tabulation of Laplace coefficients in 
dynamical astronomy, and at the time the author was unaware of any 
possible electrical application. This table appeared on reflection to offer 
a means of recalculating Spielrein’s values by a totally different method, 
namely that of numerical quadrature of the integrals occurring in (2). The 
author was in the best position to carry this out, since he possessed manu- 
script 12-decimal working values of K and E for k = 0(-01)-7(-005)1, which 
would give a more accurate basis for numerical quadrature than the 
published 10-decimal values for k = 0(-01)1. 

It was decided to tabulate J and J correct to 10 decimals, both for the 
sake of uniformity with the table (6) of K and E, and for the convenience 
of anyone who might wish to tabulate related integrals from equations such 
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as (17). In the case of f(a) there seemed no point in going beyond 6 decimals: 
this is the maximum number of decimals given (inaccurately) by Spielrein, 
and involves giving 7-8 significant figures, which are adequate for all 
electrical purposes. All quantities were to be calculated and, where 
possible, differenced to about three more figures than were to be printed. 

It was found necessary to develop series for computation near singulari- 
ties at a = 0 and « = 1, but numerical quadrature remained fundamental, 
In fact, J and J are otherwise not easily calculated in part of the range; 
the values given for them in the tables between «a = -5 and « = ‘8 rest 
entirely on quadrature, since no formulae have been evolved which enable 
isolated values of J and J to be computed to about 12 decimals much 
inside this interval, although with further labour such formulae could be 
worked out. 

Known series for small « are easily used to any number of terms. But 
the series for the larger values of « (1—a small) can be extended only with 
considerable labour. It was decided to give information regarding auxiliary 
functions in Tables 2 and 3, which relate to the case 1—« small, except that 
the last third of Table 3 relates to the case « small. 


4. Construction and description of tables 

A. Low and medium values of « 

The tables of J and J rest mainly on two highly accurate numerical 
quadratures carried out by Miss Olive E. Vanes. 

From the author’s manuscript values of K and EF, 12-decimal values of 
X K—E were formed for k = 0(-01)-7(-005)-93, and differenced. The 
value of J at « = 0 being G—}4, numerical integration (using central even 

1 

differences, up to the twelfth when necessary) gave J = | X dk for 


x 
a = 0(-01)-75(-005)-90, 14 decimals being retained. Checks were made 


at « = :1(-1):5 by the power series in «, the term in «** being included at 
a = +5, and at a = ‘88, -89, -90 by the extension of the series (9) in y 


described below. Agreement was found in all cases within 5 units of the 
fourteenth place. The 14-decimal values of J were checked by differencing. 
Table 1 gives J to 10 decimals at interval -01; differences are not printed, 
but are such that the table is interpolable by ordinary methods up to about 
« = ‘90, naturally with high-order differences near this value of «. 

In the case of J, a singularity occurs at « = 0, because the integrand 
Y = (K—E)k- has an infinite part }7k— at k = 0. Thus J+ }7lIna was 
calculated to 13 decimals for « = 0(-01)-35 from the power series in «, and 
checked by differencing; the results are given to 10 decimals at the end 
of Table 3, the corresponding values of J being inserted in Table |. 
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*12449 26518 53°178968 
O°10747 05505 55°748695 
08973 66646 58°844582 
O7105 36874 62767886 
‘O51OI 28237 68°193174 
*02873 82572 77°266993 
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TABLES OF TWO INTEGRALS 
TABLE 2 
Coefficients in Expansions of Auxiliary Functions 


[he argument in the first column indicates the power of y 














0 o’o Io 
~ ors 0°25 
4 1062 5 6406 25 0°1875 0°20312 5 
6 65 625 12044 27083 | o11718 75 0°14002 § 
8 176 80664 488 42367 0°08544 92187 10691 32487 
10 56 21338 7520 20569 0°06729 12598 008014 34937 
12 20 I< 18 24963 O°05551 52893 0°07210 57892 
14 1387 87937 5780 63778 0°04725 40855 o00199 33544 
16 56 5346 2 77040 0°04113 63691 005436 48807 
18 139 93364 1582 gs259 003642 28268 0°04540 63621 
20 54011 2 45530 0°03267 93696 0°04362 4057 
22 72353 5 78437 0°02963 42405 | 0°03970 12107 
24 )7 90755 195 47162 0°02710 86005 0°03642 53767 
| P 
I I I°o I 
2 I 2°25 9/4 
3 ac/t 108222 33323 3°22222 22222 29/9 
4 13/1 75 4°14062 5 265/64 
5 1057 25 5°04093 75 1O131/3200 
6 507/384 44791667 | 5°93315 97222 6835/1152 
= 54906 517586 6°52149 50043 
s 724 60937 7°70758 05664 
9 5Q 31522 5°59240 25935 
10 9°47939 39412 
tI . 10°35987 75585 
12 i 1124300 O7511 
So 
| = ——a 
° I O'5 1/2 
2 j 4166 66667 0°14930 55550 43/288 
4 11/288 I 94444 0°00666 66667 1/150 
6 103/10752 5 79613 090109 14802 98579/90316800 
g 25 22280 0°00028 79698 





0°00009 80036 








TABLE 


Values of Auxiliar 


x 1019], 1010] 
0°88 | 5821 88538 | 5543 69431 
"89 | 535019647 | 5117 03859 
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101° J, 


47 | 739° 93019 
346 | 6650 goo44 














0:90 4876 08247 4683 88166 5711 24858 5936 79887 
"OT | 4399 53003 | 4244 24830 | 5068 29462 | 5247 49413 
"92 | 392054990 | 3798 16289 | 4443 00742 | 4581 91515 
93 | 3439 11696 | 3345 64940 | 3834 69034 | 3939 04735 
"94 | 2955 23020 | 2886 73139 | 3242 68196 | 3317 92907 
0°95 2405 55270 2421 43204 2666 35392 2717 64832 
06 1980 06768 1949 77413 2105 10892 | 2137 33978 
‘97 | 1488 77845 | 1471 78010 | 1558 37885 | 1576 18189 
98 995 00840 987 47200 | 1025 62311 | 1033 39430 
99 498 75105 496 87151 506 32695 508 23539 
| 
1:00 ° ° ° ° 
f 52 52 
x My re) } le o* 
o-88 | 11°81439 44 +316 | 5:89900 73 —II 31 
“89 11°87689 68 311 5°93160 33 Il 13 
0°90 11°93943 03 306 5°96408 80 —1095 
‘OI 12°00199 44 3 OI 5°99646 32 10 75 
“92 12060455 56 296 6°02573 06 10 OI 
93 I2°12721 24 292 6:06089 19 10 44 
“04 12°18986 54 2 37 6°09294 88 10 28 
O95 | 12:2! 283 | 6°12490 29 —10 12 
‘96 | 12°3 278 6°15675 58 997 
‘07 12°3 274 6°18850 89g g 82 
‘98 12°4 270 6°22016 39 9 67 
‘99 12°5 2 66 6°25172 22 9g 52 
1:00 12°56637 06 +262 6:28318 53 —9g 38 
x J+47 Ina x ]+4rIna x J+41 Ino 
0:00 0°28012 83693 O-12 0°27799 81543 0-24 0°27145 d5g614 
‘OI "28011 36426 13 *27762 63862 25 *27073 90493 
-02 *28006 94570 "14 -27722 41818 26 °2 } 
‘03 *27999 57900 "15 *27079 14152 ‘27 *26913 93331 
0:04 0°27989 26319 0-16 o 0:28 0°26828 88484 
“05 *27975 99259 a7 «CI "29 -26740 41288 
06 *27959 76281 18 *30 *26648 47879 
‘07 "27940 50771 ‘19 31 *26553 04195 
0:08 | 027918 40006 0:20 | 0°27416 27855 0°32 0° 26454 05971 
‘09 27893 25145 ‘2 | *2735 “33 26351 48726 
‘IO *27865 11233 22 °272 34 -26245 27761 
‘Il *27833 97196 "23 | + # *27220 60012 35 *26135 35143 
1 
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The numerical integration started at « = -35. Values of Y were computed to 


11 decimals for a -30(-01)-70(-005)-93 and differenced. Numerical inte- 
1 

gration gave J | Y dk to 13 decimals for a -35(-01)-75(-005)-90. Series 

checks at a ‘4, -5 and at a -88, -89, -90 revealed no errors exceeding 


7 units of the thirteenth place. Table 1 gives J to 10 decimals at interval 
01; the values are interpolable by ordinary methods except at both ends 
of the range, where special provision is made. 

Values of f(a) were calculated to 9 decimals from (1) and checked by 
differencing. Values to 6 decimals are given in Table 1, together with 
modified second differences 5°, defined by 
on, 52— 0-18484-+ 0-0380868 — 0-008388 + 0-0026"... 


ind also. for the short range from « -82 to a -90, values of I4 defined 


v) 100014 — §4—0-27888-+0-0788... . 


On the use of these for interpolation, see (12); the present paper uses the 
notation I4, instead of the customary y’, in order to avoid confusion with 
the parameter y defined by (5). Everett second-difference interpolation 
may be used for 0-02 < a < 0-82, while the slight correction for [' when 
0:82 << « < 0-90 may be omitted if only 5 decimals are needed in f(a). 
The values of 82, start at a 0-02, since interpolation in f(a) for very small 
. is affected by a term proportional to o?Ina. For a < 0-02, f(x) may be 
calculated to 6 decimals from 
(1—.«)*f(a) 6-96957 043 —a3(9-0801— 30-3008 logy, «); 

» more extended expression is given in (5), p. 555, correctly except for 
trivial errors of one final unit in two coefficients (for 9-08008 and 0-33045 
read 9-08009 and 0-33046). 


B. High values of x 
For high values of «, some difficulty was at first found in calculating J 


ind J to 12 or 13 decimals. One can work out without too much difficulty 


expressions such as 


i = ces 4 i, ‘i 2) ws 
I — ~o'* 4 xy °- .in ,— x’? 1 — g/4 4 64), 
2 16 384 : x 4 64 576 
where a’? = 1—a®; such expressions, however, are not nearly powerful 


enough, unless a > 0-99. It is almost essential to use the Landenized 
rgument y defined by (5), and to extend the series (9) by use of (6) and (7). 

The power series K,, K,, E,, E, occurring in (7) are well known (7), but 
the writer is aware of no place in which the numerical values of all four 
sets of coefficients may be found. Values were therefore computed to 15 
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or more decimals, and are given to 10 decimals at the top of Table 2 as 
far as coefficients of y**, although coefficients up to those of y!° inclusive 
are all that are needed in calculating the remainder of Table 2. With the 
help of the coefficients printed it would be possible to recalculate Airey’s 
useful tables (7), which have been shown by the writer (11, p. 261) to 
suffer from last-figure errors. 

Performing on power series in y the somewhat tedious operations speci- 
fied in (6), series of the form (9), but extending to y!? and having decimal- 
ized coefficients, were worked out. The coefficients in J, and Jj, which 
apart from alternate signs are also the coefficients in J, and J,, are given 
to 10 decimals in the middle of Table 2. The coefficients up to y* check 
Spielrein’s vulgar fractions from (9), which are printed alongside; the check 
is exact, as these coefficients either terminate or exhibit simple recurrence 
within 10 decimals (15 or more decimals were used in the calculations), 
The coefficients of y’ to y!? are new. 

Values of J,, J,, J,, J, were calculated from these series to 14 decimals 
for « = -88(-01)1, and checked by differencing; the values are given to 
10 decimals at the top of Table 3, and are easily interpolable. The corre- 
sponding values of J and J obtained from (8) are given at the end of Table 1. 

Performing on power series the operations implied in (11), series of the 
form (4), but extending to y!° and having decimalized coefficients, were 
evaluated. The coefficients in S, and S, are given to 10 decimals at the 
end of Table 2, with the (amended) Lyle coefficients as far as y®, as given 
in (4), printed alongside. All these are exactly checked, although 10- 
decimal values do not demonstrate this in the case of the coefficients of »’, 
The formulae for S, and S,, given to y* by Spielrein and to y® by Lyle, are 
thus extended to y?®. 

For « -88(-01)1, it was thought best to tabulate as auxiliary functions, 
not S, and S,, but 

hi 2m(1+-a)S;, te 27(1+-«)S,, 
so that, by (3), f(«) may be computed from 


f(a) = f,In(4/y)—fe. 

Values of f, and f, were computed to 9 decimals and checked by differenc- 
ing, and are given to 7 decimals, with second differences, in the middle of 
Table 3. The corresponding values of f(«) are given at the end of Table |. 

Computations required for this paper were performed in the Mathemati- 
cal Laboratory of the University of Liverpool, a number of them by Miss 
Olive E. Vanes, B.Sc., as part of a programme approved by the Department 
of Scientific and Industrial Research. The author wishes to express gratefil 
thanks for this indispensable assistance. 
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A SIGNED BINARY MULTIPLICATION TECHNIQUE 


By ANDREW D. BOOTH 


(Birkbeck College Electronic Computer Project, 
21 Torrington Square, London, W.C.1) 


[Received 1 August 1950] 


SUMMARY 
A technique is described whereby binary numbers of either sign may be multiplied 
together by a uniform process which is independent of any foreknowledge of th; 


signs of these numbers. 


In the design of automatic computing machines it is necessary to have 
available some means of multiplying together two numbers whose signs 
are not necessarily positive. This, of course, is completely trivial when 
the process is to be performed by a human operator since a large number 
of processes exist. However, few of these seem to be suitable for mechaniza- 
tion with the types of circuit currently available on account of the com- 
plexity of the discrimination required for their execution, and the problem 
is to find a procedure. which can be engineered with the minimum of 
equipment. Several ways of accomplishing this have been used to date, 
all more or less unsatisfactory, for example: 

(1) the machine may use numbers in the form (sign) (absolute value of 
number), in which case, although multiplication (and division) are 
particularly simple, the much more frequent operation of subtraction 
needs special circuitry; 

(2) negative numbers may be represented in complementary form 
mod 2? 

when it is necessary either first to convert them to positive form, multiply, 
and then to correct the resulting product to its signed value by means of 
a special sub-routine; or to apply appropriate corrections to the product, 
obtained in the usual way, by neglecting the fact that the numbers may 
be non-positive. The nature of these corrections is seen from the following 
discussion. 

Assume that the machine in question (this is the case in the machines 

under construction in this and many other laboratories) deals with negative 
numbers by taking their complements mod 2, then: 


+m =m 


—m 2—m. 


[Quart. Journ. Mech. and Applied Math., Vol. IV Pt. 2 (1951)] 
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Whence if two numbers m and r are to be multiplied, the machine will 


generate the following results: 


L-m>x +r mr, (a) 
mx +r = 2r—mr, (b) 
mx —r 2m—mr, (c) 
mx —r = 4—2m—2r-+mr. (d) 


Thus, in order to correct (b), (c), and (d), it is necessary to apply the 
following process: 

(1) if m is negative, subtract 2r from the product obtained in the normal 

manner; 

(2) ifr is negative, subtract 2m from the product obtained in the normal 

manner 
the application of both of these corrections also gives correct results in 
case (d), since if m and r are negative, subtraction is in effect addition and 
since operations are all mod 2 the added 4 is in any case ignored by the 
machine. 

It is evident that the application of this correction process involves 
examination, by the machine, of the signs of both m and r and this, in 
turn, requires for the efficient engineering of the sequence the storage of 
the signs of m and r in auxiliary circuits (1). 

Such correction operations as envisaged above are highly undesirable, 
and it is natural to inquire whether any process exists whereby multiplica- 
tion can be performed in a uniform manner without the necessity of any 
special devices to examine the signs of the interacting numbers. That this 
was a reasonable quest was rendered probable by the development, by 
Burks, Goldstein, and von Neumann (2), of the so-called ‘non-restoring’ 
process for the division of signed binary numbers, and a somewhat com- 
plicated process for multiplication had in fact been suggested previously 
by Rey and Spencer (3). 

An extremely simple process (both mathematically and technically) has 
now been evolved and forms the subject of this note; it was suggested by the 
standard ‘shortcutting’ method of multiplication used on desk machines 
operating in decimal scale. Thus, to multiply by 057737 (i.e. +57737) on 
a desk calculating machine using the shortcutting method, the operator 
effectively multiplies by 142343. Similarly to multiply by the number 
977563 (i.e. —22437) the operator effectively multiplies by 1022443. 
The corresponding process for binary numbers gives for multiplication by 
the positive numbers 0111011, multiplication by 1001101. The process 
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is precisely this, provided the multiplication starts with the least significant 
digit, and may be described as follows: 
To multiply two numbers m and r together, examine the nth digit (m | 
of m, 
(1) Ifm, = 0, m,,, = 0, multiply the existing sum of partial products 
by 2-1, i.e. shift one place to the right. 
(2) Ifm, = 0, m,,, = 1, add r into the existing sum of partial products 
and multiply by 2-*, i.e. shift one place to the right. 
(3) Ifm, = 1, m,,,, = 0, subtract r from existing sum of partial products 
and multiply by 2-1, i.e. shift one place to right. 
(4) If m, = 1, m,,, = 1, multiply the sum of partial products by 2-1 
i.e. shift one place to the right. 
(5) Do not multiply by 2-1 at m, in the above processes. 
Note that if m is given to n digits, at the start of the process it is assumed 
that m,,,, = 0. 
In proving this process it will be assumed that operations are (mod 2). 
Thus 
mM = My.M,Mz...M,, 
My. 2°+-2-1m,+2-2mg...+2-"m,,...+2%my (m, = 0, 1). 
Now consider the multiplication of r by the number 0.0...01,, 00...0, ie. by 
2-". It is seen that the process (1)—(5) gives, for the product: 
—rxX 2-"+r x 2-"+1, 
1.€. rx 2-"(2—1) = 2-4, 
which is the correct result whatever the sign of r since the multiplication 
by 2-" is achieved by successive right-shift operations. 
Thus, up to mp», the following results will be obtained in the complete 
multiplication by m: 


(+m) xr is correct for all signs of r, 
(—m) Xr = (l—m)r = r—mr for all signs of r. 
At stage mp», however, the quantity m,).r (m, = 0, 1) is subtracted from 
these results, giving: 
(+m)xr—0xr = +r for all signs of r, 
(—m) xr—1xr = (l—m)r—r = —mr for all signs of r. 


Whence the process outlined gives the correct product for mxr (mod 2) 
whatever the signs of m and r. 

In conclusion an example of the application of the procedure to each 
possible combination of sign will render the process clear. 
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0-101 (+2) 


subtract r 
shift right 
add r 


shift right 
subtract r 


shift right 
add r 


no shift 


L-O1ll ( 5) 


subtract r 
shift right 
shift right 
add r 


shift right 
subtract r 


no shift 


0-101 (4 3) 


subtract r 
shift right 
add r 


shift right 
subtract r 


shift right 
add r 


no shift 





r 0-110 (4 


1-1010 
0-110 


0-0110 
0-001,10 
0-00110 
1-010 
1-01110 
1-101,110 


1-101,110 
0-110 


0-011,110 (mod 2) 
0-011,110 (= +433) 


32 
r = 0-110 (+ 3) 


1-010 
1-101.0 
1-110,10 
1-110,10 
0-110 


0-100.10 

0-010.010 
0-010,010 
1-010 


1-100,010 


1-100,010 (= —38) 
r = 1-010 (—3) 
0-110 
0-011,0 
0-011,0 
1-010 
1-101,0 
1-110,10 
1-110,10 
0-110 


0-100,10 

0-010,010 
0-010,010 
1-010 


1-100,010 
1-100,010 (= —}3) 





m, = 0, m, = 1 add r 0-001,10 
1-010 


1-011,10 
shift right 
m, = 1, m, = 0 subtract r 1-101,110 
0-110 


0-011,110 
no shift 


REFERENCES 
Electronic Computer (Princeton, 1947). 
Computing Instrument (Princeton, 1946). 


Conference Report, 1950. 
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(4) m = 1-011 (—8) r = 1-010 (—3) 
ms = 1 (m, = 9) subtract r 0-110 
shift right 0-011,0 
ss - i,m, = I shift right 0-001,10 


1-101,110 


0-011,110 (+38) 


1. A. D. Bootn and K. H. V. Brirren, General Considerations in the Design of ai 
2. A. Burks, H. Gotpstern, and J. VON NEUMANN, Logical Design of an Electroni 


3. T. J. Rey and R. E. Spencer, Sign Correction in Modulus Convention, Cambridy 
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THE INFLUENCE OF 
COMPRESSIBILITY IN ELASTIC-PLASTIC BENDING 
By J. W. CRAGGS (St. John’s College, Cambridge) 

[Received 2 November 1950] 


SUMMARY 
In the classical theory of bending of a uniform beam by terminal couples all stress 
components except the longitudinal one are neglected. This assumption is shown to 
e incorrect for an elastic-plastic material unless it is incompressible and the magni- 
tude of the hitherto neglected stress components is estimated by finding those arising 
in a simpler elastic-plastic problem of the same nature. 
1. Introduction 
WHEN an elastic-plastic material is stressed to the yield limit, any plastic 
flow which occurs takes place without change of volume. In many cases 
the calculation of the stresses in such a material is simplified if the material 
is taken to be entirely incompressible. Such a problem is that of the flexure 
of a beam, and it seems desirable to investigate the magnitude of the errors 


introduced by this assumption. 


2. The nature of the problem 

In his treatment of the bending of beams Nadai (1) assumes that all 
stress components other than the longitudinal tension or compression are 
zero. 

Take axes Oz along the beam, Oy as the neutral axis, and Ox in the 
direction of the principal radius of curvature. Then, for an ideal elastic- 


plastic material, the assumption is that 


ae c. Tne yy Oy, = 9, 
with C2, Exx in the elastic region, |Ex«xr| < Y, 
a Y for Exx > Y, 
and o Y for Exzx - 


where « is the curvature, E is Young’s modulus, and Y is the yield stress 
in tension and —Y in compression. 
The corresponding strain components in the elastic region aret 
( KX, Cos Cuy VKIX, Cry Cr: Cuz 0, 
and the displacements 


u 4n(vax®—vy?+-z?), v KVLY, w= ane, 
where v is Poisson’s ratio. 


The tensor definition of the shear strain components is used throughout. 


[Quart. Journ. Mech, and Applied Math., Vol. IV, Pt. 2 (1951)] 
5092.14 
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In the plastic region, symmetry and incompressibility lead to 


: 1 — 
Cnr Cyy — 2cz — ge, 
thus the displacements are 
u = —4dx(hu?—4hy?+2?)+h,, v Leay tk, w = Krztky, 


where the arbitrary functions k,, k,, k, are at most linear in 2, y, and z, 
There is thus a discrepancy in u, proportional to (4—v), which cannot be 
removed by any choice of the k’s. It follows that unless v 1. the 
incompressible case, the stress components perpendicular to Oz are non- 
vanishing and must be determined for a complete solution of the problem. 
This determination presents two main difficulties. First, the complete 
plasticity equations will contain the curvature of the beam as a third 
independent variable, and the problem reduces to a non-linear partial 
differential equation in three variables. Secondly, the position of the 
elastic-plastic boundary depends on the stress components and is thus 
unknown until the solution is found. To avoid this latter difficulty, while 
obtaining some idea of the general effect of the hitherto neglected stresses, 
Professor Sir Geoffrey Taylor suggested investigation of the following 


related problem. 


3. A simplified model 


Consider the uniform extension in the z-direction of an infinite plate, 


bounded by the surfaces y = +a and such that the material on the 
negative side of the plane 2 = 0 is perfectly elastic, while that on the 


positive side is a Prandtl—Reuss material with the same elastic constants 
but with a finite yield limit, Y. 
The problem is two-dimensional, and neglect of ,.,, ¢,,,, and a,,, leads to 


Crr e vy 


and a € i te L(g EY) (a > 0). 


ra vy 


vé (x <0) 


where € is the tensile strain in the z-direction. Thus e,,, and v, the y- 
component of displacement, are discontinuous at x = 0. It will be 
observed that in this case the hitherto neglected forces are such as to 
compress the thickness of the plate on one side of 2 = 0 and expand it 
on the other. [In the flexure problem the forces must be such as to 
produce bending in the ay-plane, of opposite signs on the two sides of the 
elastic-plastic boundary.] The discontinuities introduced are shown in 
Figs. 1(a) (bending) and 1(b) (simplified model), where the material is 


supposed elastic for x < 0 and plastic for x > 0. 
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COMPRESSIBILITY IN ELASTIC-PLASTIC BENDING 
J y 
: > x | ———_> x 
(2) (t) 
Fic. | 


4, The equations for solution 
The equilibrium equations, which hold for the whole plate, may be 


satisfied by writing 


c*x oy” Gs o*x Cx, to c*y Oxoy, (1) 
where x is a stress function independent of z. 
For 2 0, the stress o,, is then 
0,. hé t vV2x, (2) 
where V? 0° /ex* +c? cy*, € is the (uniform) strain component e,,.. Also 
, satisfies ’ 
X tist Vty 0. (3) 
For 2 0 the stress o,, is determined in terms of x by the yield con- 
dition ) > 9 9.9 ré 
ee es 302, Y2. (4) 


and the strain components are given by the flow equations 


kée..—é Cy,+ Ger) Ke ss Cay v(o,,+6--) 
+ a - ) 
ZC Oo O Oo yy Ory 0; 
hé,,—(1+v)é,, iD 5.2 VWGz7+Gyy) (5) 
» i“ ? 
0» y -0.- Ory Tyy 
where dots denote differentiation with respect to €. 
The compatibility condition. 
Vaal Fake Oe 
9 ry vy . 
: - t 5 Q, (6) 
cy” Oxcey ox* 


then gives a non-linear partial differential equation for y in terms of x, y, 
and €. The solution of such an equation offers such complications that it 
seems essential to adopt some linearization procedure. 


9. Expansion in powers of (4—v) 
It has been observed that for v } the solution reduces to x 0. It 
therefore seems reasonable to expand y in powers of 4—v x and write 
x vd +- ays +higher terms. (7) 
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The values of ¢, ys, ete., can then be found from a sequence of linear 
equations, and it may be hoped, since }—v is not large, to obtain a useful 
answer by neglect of «? and higher terms. 


Then (4) reduces to i, «= Fe, (8) 
and (5) tot 
oO, . 3a 
bv b+ 5 (Puy —bae) + Fy Puy —berx): 
; 3a |: . 3a 
Tin ao 1E (.—Pyy) + ay (¢,., Py) 
: da : Say 
and bu 7) Pry - Oy Dae (9 
Use of (6) then leads to 
¢ EK 
+—lV4%=0 (x >), (10) 
f } 
and as before Vd 0 &@< 0) (1] 
The linearized problem requires the solution of (10) and (11) with the 
] | 
conditions a om 0, ¥ +a, (12) 


the second partial derivatives 


— Pyy tend to zero as x > 0, (13) 
¢ = 0, all x,y, at € = Y/E, (14) 
and Trp Try, Cy» U, and v are continuous at 2 = 0. (15) 


Try a solution of the form 
h = F(E) P(x, y), (16) 


where F is a function of € only and ® is biharmonic in x, y and independent 
of €. This expression satisfies both (10) and (11). 
Now for x < 0 


Cuy Evé T (Sax 4)(¢,, Cunt 
to the first order in a. Continuity therefore demands that 


3 eo aod 
LB : 
4K | . 0 


Ox? Oy? 


‘ D ‘ 9 
oo oF 


BF (EPO CD _ BFE) [eP@ - = 
4h : JIr=0 ar = Jz=0 


where accents denote derivatives. 


<a 
an" Oy 


+ The last fraction in (5) is E/2Y+O(a?). 
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near 20 2p ” 
efi This implies 5 . const = A, say, (17) 
, On" ad 
CeDp ao 
and : const B, say. (18) 
ot a OY 1 x=0 
Then (B—A)F’(é)+(E£/Y)BF(€)—(4E/3) = 0 
. 1 (EB(é—Y/E)) 
hence F(€) ] “xp a (19) 
_— = Bl P\-Y(A—B) | 
Now (17). (18), and the boundary conditions (12) are satisfied by 
@ Ax?/2+¢6* (4% < 0) 
. Bau?/2+4* iz > ¥), 
where 4* is a solution of V4d6* 0 with 
*, 0. d* A aty Lae, « < 0, (20) 
() d* ’ Oo. Ad* B at y +-a, x 0. (21) 
1] * may now be determined by integral transform methods ast 
1+ B 
h the + ™ 
2 
(12 
B—A f[ AysinhAysinhAa—cosh Ay(Aa cosh Aa+-sinh Aa) _ sinAx dA 22) 
. 2Aa+-sinh 2Aa r3 i 
(13 0 
and the conditions at infinity and the condition ¢*, —¢*, > 0 as x > # are 
14 


satisfied when A+ B 0. Thus 


(15 )} | ERE s a 2 | , 
f l—exp = + I (23) 
| <P(5 sah 2 2 | 


r where the upper sign is taken for x > 0 and the lower for x < 0, and J is 
the integral in (22). 
ndent } 
| 6. Results 
The stress components have been computed for 0 < y/a <1 and 
ra<2. Tables 1 to 4 give quantities s,, s,, s,, and ¢ such that 
| i Z uv Z 
o,,/Y (4—v){1—exp($— HE/2Y)}s, 


with similar relations between s,, t and a,,,, Cam and 


o,,/Y l+(4—v){l—exp(4— HE/2Y)\s, (x > 0) 


BE/Y —2v(4—v){1—exp(}— EE/2Y)}s, (a < 0). 


This integral is divergent at the lower limit, but the integrals obtained from it by double 
ferentiation under the integral sign are convergent. (22) may therefore be regarded as a 
mnemonic for the 


stresses. 
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0°00 
0°25 
0°50 
O75 | 
I'oo | 


, 8, > 0-667. 


TABLE 4. tx 103 








where at x 


7. Conclusions 

i. The mean value of s, over the table is 0-031. The total force required 
for the extension of the composite specimen is therefore given sufficiently 
accurately by neglect of the forces in the cross-section plane. 

ii. The stresses in the plane rise to the order of (4- 
y =a. The area over which such values are reached is small and it seems 
entirely. 
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reasonable therefore in this problem to neglect a,.,, 
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However, in the bending of an elastic beam these components may have 
a considerable effect on the position of the elastic-plastic boundary. It is 
hoped that a further investigation at present in progress may throw more 


light on this question. 
I iii. The method used here can be easily extended to the extension of 
composite prismatic beams of other cross-section, but in view of the 
artificial nature of the problem the work is scarcely worth while. 

In conclusion I would like to express my thanks to Professor Sir Geoffrey 
Taylor for his constant help and encouragement in this work. 
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SUMMARY 

Following the methods outlined in a previous paper, complete solutions are given 
for the stresses in an aeolotropic elliptic plate under the most general type of edge 
loading in the plane of the plate. Since the branch points of the transformations used 
lie inside the boundary of the plate care must be taken in choosing the potential 
functions to ensure that these are single-valued. The general method of dealing with 
all such cases is illustrated. 

1. Introduction 

In the previous paper of this series by Livens and Morris (1)—subse- 
quently referred to as (1)—complete solutions have been given to the 
problems relating to (i) the effect of an elliptic hole on an otherwise uniform 
distribution of stress in an infinite aeolotropic plate, and (ii) the stress 
distribution produced in an infinite plate by specified forces at a point or 
on the edge of an elliptic hole. ° 

It was pointed out in (1) that the particular case of the elliptic hole in 
the infinite plate does not introduce the special difficulties associated with 
the branch points of the transformations used, which lie outside the 
boundary of the hole, and which do arise in the consideration of problems 
connected with holes of more general shape in an infinite plate. 

The special difficulty of these problems arises from the fact that the 
potential functions involved might, unless carefully chosen, change in 
value in crossing any cut joining two such branch points of the transforma- 
tion used. These difficulties, however, are not confined to problems of 
holes in an infinite plate, but do in fact arise when we consider the alterna- 
tive problem of determining the stress distribution in a finite plate whose 
edges are subject to specified loading in the plane of the plate. In such 
cases we are not concerned with the infinitely-distant parts of the plane. 
so that the convergence of the potentials there is not a condition of the 
problem. We have, however, to determine a potential function Q which 
is not only regular over the area of the plate, but also determines regular 
stresses, and satisfies a specified boundary condition. 

Thus the transformation 


2= c cos(€+i«), C = €+77, (1.1) 


[Quart. Journ. Mech. and Applied Math., Vol. IV, Pt. 2(1951)] 
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is such that the boundary corresponding to 7 = 0 is an ellipse, but when 
this transformation is used for the finite elliptic plate, we find that the 
usual transformations arising in the analysis are such that the branch 
points lie within the plate. The particular difficulties mentioned above 
will therefore arise. The problem of the finite elliptic plate under specified 
edge loading is in fact one of the simplest cases of this type, and it is 
therefore proposed to indicate the approach to such problems by examining 
three cases of typical edge loading for the elliptic plate. 

Solutions of this problem when the material of the plate is isotropic have 
| been given by Scherman (2) using the general methods of Muschelisvili and 
Kolosoff. Stevenson (3) also gives methods which can be applied to the 
problem of the finite isotropic elliptic plate. 


2. The stress potentials 


Following the methods of (1) we have to determine a stress function 





which takes the form 


Sloe log 9 
Q2 > Se m(Zm) + ¥ Sen (=m) f> (2.1) 
i | m 
where, as usual, z,, = z-+A,, Z. 
Using the transformation z = ccos({+ia) we have 
Zm == €C08(C+ia)+A,, Ecos(C—ia), (2.2) 
and then on the boundary of the plate, where 7 = 0, we have 
m= %plt) = $(ce*+Ay e-*)(t-+-y23,/t), (2.3) 
where ¢ = exp(—ié) and 
yi, = (ce-*+Ay, Ge%)/(ce* Ay, Ce-*). (2.4) 


In general, this expression determines ¢ as a function of z,, which reduces 


to exp(—ié) on the edge of the plate. It is a transformation from the 


| z,,-plane to the t-plane such that the distorted ellipse in the z,,-plane 
corresponding to the boundary of the plate in the z-plane is transformed 
to the circle |f | in the ¢-plane. Moreover, we have 

| z(t) k(ce*-+A,, €e-*)(1—y?, /t?), (2.5) 
and this will be zero at the points t = +y,,. These two points are the 

| branch points of the transformation and correspond to the foci of the family 
of ellipses in the z,,-plane. If A,, < 1 so that |y,,| < 1 these points corre- 
spond to points inside the circle |t 1, and in fact the degenerate ellipse 


of the family—the double segment joining the two foci—corresponds to 
the circle of radius |y,,| in the t-plane. 


These two points being the branch points of the transformation we have 


now to choose the corresponding part Q,,,(z,,) of Q to be such that its value 
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is unchanged in crossing any cut joining them. We secure this by noting 
that for the transformation 


2m = Z(t) : 3(ce*+A,, Ce-*)(t+-y7, t), (2.6) 


m 
there are two values of ¢ for every value of z,,. If t = t,, be one of these, 
we can write 


em = 3(ce es i Am ce-*)(t-+ Yin t) = $(ce*+-A,, €e-*)(tm+Yin tm), (2 


“m 
so that the other root is y?,/t,,. 


Taking a general case, suppose that a characteristic boundary condition 
for the potential function Q is of the form 


y Bt. (2.8) 
r x 


Then we have already seen in (1) that a suitable form for Q,,(t,,) is 


+ CO 


r 9 
3 Amt ? (2.9) 
x 


r 


where the a,,,,. are constants to be determined from the boundary conditions. 


mr 
But for any function of t,, we know that in going round either of the two 


branch points, the two roots ¢,, and y?,/t,, are interchanged, and therefore 


7 


by writing $a 
On (tm) a > Ome (tin +Yin tn)» (2.10) 


; 
we can ensure that in such a passage round one of the branch points, the 
value of Q,,(¢,,) remains unchanged. 

This is then the general result and the remainder of this paper is devoted 
to its application to problems in which the boundary value of Q takes 
particular forms, including one in which the stresses at the edge of the 
plate are non-equilibrating. 

3. The simplest type of edge loading is that in which the stresses are 
constant and given by 

qn—ing = —(P+iQ) = —Dp, say. (3.1) 
These stresses on the edge of the disk, however, are equivalent to a force 
at the origin having components X, Y given by 


X+tY =1 [ (nn—iné) dz, (3.2) 
C 
and a couple about the origin of moment 
G = re [ (nn —iné)z dz, (3.3) 
o 


the integrals being taken around C, the boundary of the disk. With the 
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above value for 7n—iné we notice that X, Y, G will all be zero only if p, 
is real, and, as the simplest example, we shall therefore assume 


1) 7) r. né 0. 


il 


The boundary condition (cf. (1), section 6) 


é 
o=2( (m-it)@ a (3.4) 
ad és 17) - VT 9 v. 
; 17) Is dé” 
then becomes Q Pc(e*t+e-t-1), t = e-*, (3.5) 


An appropriate form for Q is therefore 
ppro} 


Diao Bite 7” 
Q, (t,)+Q.(t 2) +7 Mali) + 7 Qaltr), (3.6) 
1 2 
W here Qinltn) (ce "| An C *)(Ay_t m +73, Am = ), (3.7) 
the constant (ce*+-A,,ée-*) being inserted for analytical convenience. On 
the boundary, t,, = t;,1 = t and therefore, —_ 
2 ce t+, 
7m cet A, Ce —a? 


the boundary condition is satisfied if 


- ‘ a, ,. de ,- 
a,(ce*+A, Ge-*)+-a,(ce* +A, Ge-*) + x, (ee “A, ce*) + 57 (6e-*+ A, ce*) 
1 2 
—Pce*, (3.8) 
and 
a.(ce-*-+. Fe e-% 1). Fen) ay Fe*+-J. ce-%) 4 dy a LA, a 
1(Cé tA, Ce )+- a,(ce T Ag Cé +5 Ay )+ 5 (eet on ce ) 
1 2 


= —Pce-*, (3.9) 
These two equations, together with the two conjugate equations, determine 
a, and dg, viz. 
l Pr,(1+A3 l PA,(1 
a, a a a ode (3.10) 
2 (Ay —A)(1—A, Ap) 220, I-AA) 
As usual we find the hoop stress [£é é], 9 at ol edge of the plate by 
evaluating [€£+-77],_», that is 


which in this case is easily seen to be —2P. Thus, since [77],-0 = —P, 
the hoop stress is [éé =0 r. 

The hoop stress at the edge of the plate is therefore independent of the 
aeolotropy of the material of the plate, and is also the same as for a circular 
plate. This is true independently of the orientation of the principal axes 
of the ellipse relative to the principal axes of stress of the material. 
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4. Taking as a second example the simplest case of a variable loading 
of the edge defined by 
mn—ink = pye~+p_e*, (4.1) 
we find that the boundary condition is 
Q C(p,e*—p ,€~“)log t | k(ce “nt 21 ce¢ “Dp  ). (4,2) 
The presence of the logarithmic term in the boundary condition for Q 
arises from the fact that the forces applied to the edge do not form a 
system in equilibrium, and are in fact equivalent to a single force at the 
origin having components X, Y given by 


r .>, dz 


X+iY i | (nn—in€é) = 
dé 


d& = cn(p,e*—p_,e-*). (4.3) 
0 
To maintain the plate in equilibrium it would be necessary to apply a 
balancing force at some point in the plate, and in this case, therefore, we 
can only find a solution to this modified form of the problem. 
We have seen in (1) that the complex stress potential due to an isolated 
force P applied at any point z, in the plate, in a direction making an angle « 
with the 2-axis, is given by 


4 
p 3 Ae log (2 -2 mo): (4.4) 


m m=1 


where Z,,9 = %+A,, 2%, and the constants A,, are given by 


7 


, ; Ss S49— Gua s 
A, he _ {812 ne O22 cos € { pics 22 i sin ¢| 
(a, X»)Syo 77 1—A, 1+, 
A . A, P . [S12 M1922 ogg ¢ 1-12 Ye P22 sine (4.5) 
: (41a )Sg97{ 1 2 1+-Ay 
A, A, r,; A, A,/d, 
A As Pe‘ : 
and are such that A+ hy SF x . (4.6) 


7 A, As 7 
In the above case it is obvious that the stresses at the edge can be 
balanced by a force at the origin, and we can take z, = 0. When this is 
so, and writing again 


oa ace ie - 1.2 o 

Zm\tm) - 3(ce*-+A,,, 6 *) (tin TYm bm) (4.7) 

ee > _— % a . 
we have Qom = A m log(t,, “You t..) —+-eonst. (4. ) 
In the neighbourhood of the edge of the disk where |t,,| = 1, since |7| <!. 


this can be expanded in the form 


— — +1 
A), > = y2nt "+ const., (4.9) 
n 


Q 


= A,, logt 


Om mT 


and 


0, (¢ 


m 
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and 
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and we now add the appropriate terms to form the corresponding function 


Q (t,,) required to satisfy the boundary condition. 


m 


Remembering the general result we write 


x 


“A "4 ‘ ‘ 
Qn (tn) = Qom + > (—)"2* (B+ yt), (4.10) 
=? 
ind the condition on the boundary will be satisfied if 
i, Ae , : : 
A, 1, c(pyé S20", (4.11) 
1 As 
Cir -4 , Cor 4 A, -» A, —2 | In x 9 
. Cg Y2 Y1 Y2T seP-1€ (4.12) 
11 1 Ay l \, A, / As / 1 
Cy, Con | x wt = ot 2 oe x ‘ 
\ yt un Yat a1 Ye A, YitAsyet+aepie™. (4.13) 
l | 
und if. for x | P 
( C 4 { 
l —4) 2n -4n “*1 an #*2 —an 
Cin te ; yA pin 4 2 pan (4.14) 
l 2 l As 72 A, 71 As 2 
Cy C2, = s4nj= An {_ 2n {2m 4.15 
CinY1 Con V2 1V1 2¥2 (4.15) 
i i 
Equations (4.6) and (4.11) give immediately 
Pet mc(p, e*— p_,e-*) (X+iY), (4.16) 


which simply expresses the result that the force at the origin must balance 
the resultant of the stresses round the edge. The other equations, together 
with the conjugate equations, are sufficient to determine the coefficients 
| Cin: Cg, for all the values of n 
| We denote by A, the determinant of the coefficients, namely 


2r r 
1 


m7 
por 


~~ 
! 


| ; 


> 
> 


tw 


oe 


4 | 
rr 


vr yey 

. = 
. = \ery(Ay—Ag\* wee ~ \ory (1 —A, Ae\? 
L+ (9171 ¥2 2)” ( XA ‘] (4171) (ra Fa)” XA, — 


Uy Ya)" + (Vr ¥2)"7}(1 —AZV(I —A3) AiAS. (4.17) 
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and by A,,, Aj», A,s, A, the cofactors of the elements of the first column 


of A,, i.e. 
(A;—A,) , ,- >, (1 —A3) ~ \op(1—A, Av) | 
A - : =~ * ad foe a a 
r NB (V1 Y2) % (272) | 
‘ | 
(A, —A,) : 9» (1 —A3) ~ \op (1 —A, Ao) 
A > 2 = ms > aad = ‘i » Yo ; 
r2 A; A. (3 1Y2) A, 2 (Yo 2) i dz | 
Ayo), <9, (123) IE 
~o, (1 ‘7 5 aT 2 | 
A. —y See of ee : Ct 
r3 1 r, A> 72 2 (Y¥1 ¥2 2) A, 42 
-o.(1—A,Ao) -., (1 —A3) . ~ \opAy—Aso 
A. wae 1°*2 _4 2r oe’ ; Vo or°*1 - 
r4 71 , 2 2 A, 22 (V1 ¥2V2) A, As 


The values of the coefficients are then given by 


A. 
o ve +4p, e-*)Ao.-4 ( y2 4 2521165 _,e*]A,,, (4.19 
72 CP ) 23 i 1, ZC] 1 24 ( } 
and, for n ~ 1, 
: ee ee 7 x o 
“in = i" + *78") ans + (Ay 7" + Agye")Ay 


‘ 9 9 A 9 . A 9 
H(A, yi"+Agy")Aons + (Ste 3 - 
Ay 2 


with A,, y, replaced by Ag, y. throughout, for ¢y5.¢,,,. 


The hoop stress in this case is found from the result that (€£+ 77) is 


re4 s Fa | 5 ty? )(cex+A ce-~) 4 


m 


mom m 
Mt 1 
9 . 4n 4n\ /2n—1/ #2 ase “pa _} np —% 9 
> =( )"Cmn (bn Ym ) bn (tn Ym )(ce Xn )}. (4.21) 
n=1 
Notice that it remains finite at the branch points ¢,, LY: 


5. The only other case that involves an analysis similar to the above 
is the one in which the resultant stress is zero, but the stresses give rise 
to a resultant couple. Such is the case when the stresses on the edge of 
the disk are given by 


mn—ink = p,e*—p_.e-Ré, (5.1) 


when G is re }c?i(p,+p_,) and the boundary condition takes the form 
Q.,=0 ce-“p_»t—ce%pyt-1+-4ce“p_, t?+-4ce-“pot-. (5.2) 
Here there is no logarithmic term in the boundary condition for Q, but 
it is obvious that the resultant couple G must be balanced either by equal 
and opposite forces P at two points in the plate, or by an isolated couple 
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olumn | of moment —G, at the origin, applied physically perhaps by a clamp. The 
potential function Q, for such a couple at the origin is easily seen to be 


Q2, dit fo & a . — . \ (5.3) 
Qniz,' 2% 2% 2) 


the stresses round a small circular hole centre the origin being equivalent 
4.18 to a couple of moment —G, with no resultant force. In addition Q, is 
non-cyclic so that the corresponding displacement w, would also be non- 
eyclic. 
We then have as usual 
Q iGX,,,/{2(ce*+-A,, Ce-*) (ty, +y>,/tm)}: (5.4) 


Om mn m 





and in the neighbourhood of the edge of the plate this can be expanded 


in the form , 
(1GAyq|77(ce*+-Ay Ge-8)} SY (—)Myzne-2n+0, (5.5) 
n=0 
1] Proceeding as before to write 
4 lp 2n+1 n,AN+2 /p2n41 re 
Qn om > )"Cmnktm TYm bn ), (9.6) 
n=0 
we choose the coefficients c,,, so that the boundary condition is satisfied. 
The above analysis disposes, therefore, of all cases in which the stresses 
give rise to a resultant force and a resultant couple. 
(4.2 
6. A more general example is one in which the edge stresses are defined 
by = ¢ ° 
7 | nn—ing = p,e™=+p_,e-mé, (6.1) 
} » 
where we shall assume that » + 1, 2 so that the stresses are self-equili- 
brating. 
The boundary condition for Q is now 
(4.2 
Q,.9 = —ic | (p,e"€—p_, e-mé (er —e-a4#€) de 
re ‘D ,t" l re “py, t n—1) ce “p , t+ CE Dp, t (n+1) (6 2) 
adove i: + a g Ps 
4 n—!] n—1 n+] n+1 
ive ris 
on and the appropriate form for Q,,(t,,) is therefore 
ajuin—1) 2(n +1) 
Qing Lyn—1_, Ym Am, Ly—(n—1) | Am,n +1 yn+1 4 Ym — Fmn+1 4- (n+1) 
5.1 n l n ] ‘i n+] - n+] = 
: (6.3) 
orm | The boundary condition is satisfied if 
5.2 —2r7a —2r 7 
yer yer 
a,,+a,,+ iru, Veter _ RP (6.4) 
rQ, but , A; Az 
yy equa a a 
y Lr 2r | a lrn = 27 = ’ 7 
} Ay p+ V9" S (6.5) 


| coup 
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where r = n—1, n+-1, and 


ees —ce-“p_»» Risa = ce“p_» | 


(6.6) 
S,-1 —ce“p,, S43 = ce-“p,, 

These four equations together with conjugate equations are sufficient to 
determine the coefficients @, ,, 41, %,,—-1> G41» 42,n-1 and their conjugates, 
They can be expressed in terms of the determinant and its cofactors, 


defined in section 4, as follows 
n—-1,1 T Ce-“D,, A, 1,3 +ce ~ A, -1,3 7 ce Pans An-w 
(6.7 


x =» OL _ pet __ fe-an 
AnH — Ce Pp n Aasia—& Dy Ansis cE Pn Anite ce B-n Anse 


(6.8) 


np a 
Ain 1 ce"p i 


with A,, y, replaced by A,, y. throughout for a,,,_, and ag, ,,. 
The hoop stress is again deduced from the result that €€-+- 77 is 


f—) 


9 
S n§z2(n+1)__, ,2(M4-1)0 _J -(m —2)642(m—1) _ 2(n—1)) 

. “iy = +1 bm vem —YVm j Amn —1 tin ven Ym | 

re 4 — Pd 

m=1 


(ce*-X,, Ce-*)(t2, —y?,) 


m 
(6.9) 
By inserting a sign of summation and including terms of the type dis- 
cussed in sections 4, 5 we can generalize this result to give that for the 
most general type of loading at the edge, with a balancing force at the 
centre, or a balancing couple, if necessary. 
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